From 9ac8bb2d27f01b6d9dcaa964c5662c7ec579e42c Mon Sep 17 00:00:00 2001 From: jodavies Date: Tue, 18 Nov 2025 12:11:57 +0000 Subject: [PATCH 1/3] feat(grcc): revert diawrap.cc changes from 2aba64bb56a44cb8f154e42c3d2ac779e26d8fe3 2aba64bb56a44cb8f154e42c3d2ac779e26d8fe3 itself reverted diawrap.cc to an older version. Undo those changes. Additionally remove all diagrams tests, these will be replaced later. But for now we need a commit which passes the test suite. --- check/examples.frm | 125 -------------- sources/diawrap.cc | 407 +++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 371 insertions(+), 161 deletions(-) diff --git a/check/examples.frm b/check/examples.frm index 4162b757..b1ac873a 100644 --- a/check/examples.frm +++ b/check/examples.frm @@ -1830,129 +1830,4 @@ Print; assert result("aPLUSbTO2") =~ expr("b^2 + 2*a*b + a^2") assert result("aPLUSbTO3") =~ expr("b^3 + 3*a*b^2 + 3*a^2*b + a^3") *--#] ExtComm_1 : -*--#[ Diagrams_1 : - Vectors Q1,Q2,p1,...,p8; - Set QQ:Q1,Q2; - Set PP:p1,...,p8; - #define LOOPS "2" - Local F = topologies_(`LOOPS',2,{3,},QQ,PP); - Print +f +s; - .end -# TODO: enable it -#pend_if true - assert succeeded? - assert result("F") =~ expr(" - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,p1,-p3)* - node_(4,p2,-p4,-p5)*node_(5,p3,p4,p5) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,p1,-p3,-p4)* - node_(4,p2,p3,-p5)*node_(5,Q2,p4,p5) - ") -*--#] Diagrams_1 : -*--#[ Diagrams_2 : - Vectors Q1,Q2,p1,...,p8; - Set QQ:Q1,Q2; - Set PP:p1,...,p8; - #define LOOPS "2" - Local F = topologies_(`LOOPS',2,{3,4},QQ,PP); - Print +f +s; - .end -# TODO: enable it -#pend_if true - assert succeeded? - assert result("F") =~ expr(" - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,Q2,-p1,-p2)*node_(3,p1,p2,-p3,p3 - ) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,p1,-p3)* - node_(4,p2,p3,-p4,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,p1,-p3)* - node_(4,p2,-p4,-p5)*node_(5,p3,p4,p5) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,-p3,-p4)* - node_(4,p1,p2,p3,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,p1,-p3,-p4)* - node_(4,Q2,p2,p3,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,p1,-p3,-p4)* - node_(4,p2,p3,-p5)*node_(5,Q2,p4,p5) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2,-p3)*node_(3,Q2,p1,p2,p3 - ) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,-p1,-p2,-p3)*node_(3,Q2,p1,-p4)* - node_(4,Q1,p2,p3,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,-p1,-p2,-p3)*node_(3,p1,p2,-p4)* - node_(4,Q1,Q2,p3,p4) - ") -*--#] Diagrams_2 : -*--#[ Diagrams_3 : - Vectors Q1,Q2,p1,...,p8; - Set QQ:Q1,Q2; - Set PP:p1,...,p8; - #define LOOPS "2" - Local F = topologies_(`LOOPS',-2,{3,4},QQ,PP); - Print +f +s; - .end -# TODO: enable it -#pend_if true - assert succeeded? - assert result("F") =~ expr(" - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,Q2,-p1,-p2)*node_(3,p1,p2,-p3,p3 - ) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,p1,-p3)* - node_(4,p2,p3,-p4,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,p1,-p3)* - node_(4,p2,-p4,-p5)*node_(5,p3,p4,p5) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,-p3,-p4)* - node_(4,p1,p2,p3,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,Q2,-p3,-p4)* - node_(4,p1,p3,-p5)*node_(5,p2,p4,p5) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2)*node_(3,p1,-p3,-p4)* - node_(4,Q2,p2,p3,p4) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,Q1,-p1,-p2,-p3)*node_(3,Q2,p1,p2,p3 - ) - + node_(0,-Q1)*node_(1,-Q2)*node_(2,-p1,-p2,-p3)*node_(3,p1,p2,-p4)* - node_(4,Q1,Q2,p3,p4) - ") -*--#] Diagrams_3 : -*--#[ Diagrams_4 : - Vectors Q1,Q2,p1,...,p17; - Set QQ:Q1,Q2; - Set PP:p1,...,p17; - #define LOOPS "6" - Local F = topologies_(`LOOPS',-2,{3,},QQ,PP); - .end -# TODO: enable it -#pend_if true - #time_dilation 2.0 - assert succeeded? - assert nterms("F") == 2793 -*--#] Diagrams_4 : -*--#[ Diagrams_5 : - #define LOOPS "6" - Model PHI3; - Particle phi,+1; - Vertex phi,phi,phi:g; - EndModel; - Vector Q1,Q2,p1,...,p{3*`LOOPS'}; - Set QQ:Q1,Q2; - Set pp:p1,...,p{3*`LOOPS'}; - Set empty:; - Local F = diagrams_(PHI3,{phi,phi},empty,QQ,pp,`LOOPS', - `OnePI_'+`NoTadpoles_'+`Symmetrize_'+`TopologiesOnly_'); - .end - assert succeeded? - assert nterms("F") == 2793 -*--#] Diagrams_5 : -*--#[ Diagrams_6 : -* #define LOOPS "6" -* Model PHI3; -* Particle phi,+1; -* Vertex phi,phi,phi:g; -* EndModel; -* Vector Q1,Q2,p1,...,p{3*`LOOPS'}; -* Set QQ:Q1,Q2; -* Set pp:p1,...,p{3*`LOOPS'}; -* Set empty:; -* Local F = diagrams_(PHI3,{phi,phi},empty,QQ,pp,`LOOPS', -* `OnePI_'+`NoTadpoles_'+`TopologiesOnly_'); -* .end -* assert succeeded? -* assert nterms("F") == 4999 -*--#] Diagrams_6 : diff --git a/sources/diawrap.cc b/sources/diawrap.cc index 723acc9c..c6d4ddb4 100644 --- a/sources/diawrap.cc +++ b/sources/diawrap.cc @@ -14,7 +14,7 @@ typedef struct ToPoTyPe { WORD *vertmax; Options *opt; int cldeg[MAXPOINTS], clnum[MAXPOINTS], clext[MAXPOINTS]; - int cmind[MAXPOINTS], cmaxd[MAXPOINTS]; + int cmind[MAXLEGS+1],cmaxd[MAXLEGS+1]; int ncl, nvert; } TOPOTYPE; @@ -199,6 +199,21 @@ int ReConvertParticle(Model *model,int grccnum) } // #] ReConvertParticle : +// #[ numParticle : + +int numParticle(MODEL *m,WORD n) +{ + int i; + for ( i = 0; i < m->nparticles; i++ ) { + if ( m->vertices[i]->particles[0].number == n ) return(i); + if ( m->vertices[i]->particles[1].number == n ) return(i); + } + MesPrint("numParticle: particle %d not found in model",n); + Terminate(-1); + return(-1); +} + +// #] numParticle : // #[ ProcessDiagram : void ProcessDiagram(EGraph *eg, void *ti) @@ -230,7 +245,8 @@ void ProcessDiagram(EGraph *eg, void *ti) // // Now get the nodes // - for ( i = 0; i < eg->nNodes; i++ ) { + if ( ( info->flags & NONODES ) == 0 ) { + for ( i = 0; i < eg->nNodes; i++ ) { // // node_(number,coupling,particle_1(momentum_1),...,particle_n(momentum_n)) // @@ -245,7 +261,7 @@ void ProcessDiagram(EGraph *eg, void *ti) // function for when we want to work with counterterms. // if ( !eg->isExternal(i) ) { - afill = fill; *fill++ = 0; *fill++ = 0; FILLARG(fill) + afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill) cfill = fill; *fill++ = 0; intr = eg->nodes[i]->intrct; for ( j = 0; j < model->interacts[intr]->nclist; j++ ) { @@ -273,8 +289,10 @@ void ProcessDiagram(EGraph *eg, void *ti) edge = eg->nodes[i]->edges[j]; vect = ABS(edge)-1; *fill++ = 6+FUNHEAD; - - *fill++ = ReConvertParticle(model,eg->edges[vect]->ptcl); // code of particle_j + int a; + if ( edge < 0 ) { a = ReConvertParticle(model,-eg->edges[vect]->ptcl); } + else { a = ReConvertParticle(model,eg->edges[vect]->ptcl); } + *fill++ = a; *fill++ = FUNHEAD+2; FILLFUN(fill) *fill++ = edge < 0 ? -MINVECTOR: -VECTOR; @@ -287,26 +305,131 @@ void ProcessDiagram(EGraph *eg, void *ti) *fill++ = 1; *fill++ = 1; *fill++ = 3; } startfill[1] = fill-startfill; + } + } + if ( ( info->flags & WITHEDGES ) == WITHEDGES ) { + for ( i = 0; i < eg->nEdges; i++ ) { + int n1 = eg->edges[i]->nodes[0]; + int n2 = eg->edges[i]->nodes[1]; +// int l1 = eg->edges[i]->nlegs[0]; +// int l2 = eg->edges[i]->nlegs[1]; + + startfill = fill; + *fill++ = EDGE; + *fill++ = 0; + FILLFUN(fill) + *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge +// + *fill++ = ARGHEAD+FUNHEAD+6; + *fill++ = 0; + FILLARG(fill) + *fill++ = 6+FUNHEAD; + int a = ReConvertParticle(model,eg->edges[i]->ptcl); + *fill++ = a; + *fill++ = FUNHEAD+2; + FILLFUN(fill) + *fill++ = -VECTOR; + // Look up in set of internal momenta set + if ( i < info->numextern ) { // Look up in set of external momenta + *fill++ = SetElements[Sets[info->externalset].first+i]; + } + else { // Look up in set of internal momenta set + *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; + } + *fill++ = 1; *fill++ = 1; *fill++ = 3; +// + *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from + *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to + startfill[1] = fill - startfill; + } + } + if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) { + for ( i = 0; i < eg->econn->nblocks; i++ ) { + startfill = fill; + *fill++ = BLOCK; + *fill++ = 0; + FILLFUN(fill); + *fill++ = -SNUMBER; + *fill++ = i+1; + *fill++ = -SNUMBER; + *fill++ = eg->econn->blocks[i].loop; +// +// Now we have to make a list of all nodes inside this block +// + int bnodes[GRCC_MAXNODES], k; + WORD *argfill = fill, *funfill; + *fill++ = 0; *fill++ = 0; FILLARG(fill) + *fill++ = 0; + for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0; + for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) { + bnodes[eg->econn->blocks[i].edges[k][0]] = 1; + bnodes[eg->econn->blocks[i].edges[k][1]] = 1; + } + for ( k = 0; k < GRCC_MAXNODES; k++ ) { + if ( bnodes[k] == 0 ) continue; +// +// Now we put the node inside this argument. +// + funfill = fill; + *fill++ = NODEFUNCTION; + *fill++ = 0; + FILLFUN(fill) + *fill++ = -SNUMBER; *fill++ = k+1; + numlegs = eg->nodes[k]->deg; + for ( j = 0; j < numlegs; j++ ) { + edge = eg->nodes[k]->edges[j]; + vect = ABS(edge)-1; + *fill++ = -VECTOR; + if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta + *fill++ = SetElements[Sets[info->externalset].first+vect]; + } + else { // Look up in set of internal momenta set + *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; + } + } + funfill[1] = fill-funfill; + } + *fill++ = 1; *fill++ = 1; *fill++ = 3; + *argfill = fill - argfill; + argfill[ARGHEAD] = argfill[0] - ARGHEAD; + startfill[1] = fill-startfill; + } + } + if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { + for ( i = 0; i < eg->econn->nopic; i++ ) { + startfill = fill; + *fill++ = ONEPI; + *fill++ = 0; + FILLFUN(fill); + *fill++ = -SNUMBER; + *fill++ = i+1; + *fill++ = -ONEPI; + for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) { + *fill++ = -SNUMBER; + *fill++ = eg->econn->opics[i].nodes[j]+1; + } + startfill[1] = fill-startfill; + } } // // Topology counter. We have exaggerated a bit with the eye on the far future. // - if ( info->numtopo < MAXPOSITIVE ) { + if ( info->numtopo-1 < MAXPOSITIVE ) { *fill++ = TOPO; *fill++ = FUNHEAD+2; FILLFUN(fill) - *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo); + *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo-1); } - else if ( info->numtopo < FULLMAX-1 ) { + else if ( info->numtopo-1 < FULLMAX-1 ) { *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+4; FILLFUN(fill) *fill++ = ARGHEAD+4; *fill++ = 0; FILLARG(fill) *fill++ = 4; - *fill++ = (WORD)(info->numtopo & WORDMASK); + *fill++ = (WORD)((info->numtopo-1) & WORDMASK); *fill++ = 1; *fill++ = 3; } else { // for now: science fiction *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+6; FILLFUN(fill) *fill++ = ARGHEAD+6; *fill++ = 0; FILLARG(fill) - *fill++ = 6; *fill++ = (WORD)(info->numtopo >> BITSINWORD); - *fill++ = (WORD)(info->numtopo & WORDMASK); + *fill++ = 6; *fill++ = (WORD)((info->numtopo-1) >> BITSINWORD); + *fill++ = (WORD)((info->numtopo-1) & WORDMASK); *fill++ = 0; *fill++ = 1; *fill++ = 5; } // @@ -329,19 +452,54 @@ void ProcessDiagram(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; +// MesPrint("<> %a",newterm[0],newterm); + Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; } // #] ProcessDiagram : +// #[ fendMG : + +Bool fendMG(EGraph *eg, void *ti) +{ + DUMMYUSE(eg); + DUMMYUSE(ti); + return True; +} + +// #] fendMG : // #[ ProcessTopology : Bool ProcessTopology(EGraph *eg, void *ti) { // // This routine is called for each new topology. +// New convention: return True; generate the corresponding diagrams if needed +// return False; skip diagram generation (when asked for). // TERMINFO *info = (TERMINFO *)ti; +#ifdef WITHEARLYVETO + if ( ( ( info->flags & CHECKEXTERN ) == CHECKEXTERN ) && info->currentMODEL != NULL ) { + int i, j; + int numlegs, vect, edge; + for ( i = 0; i < eg->nNodes; i++ ) { + if ( eg->isExternal(i) ) continue; + numlegs = eg->nodes[i]->deg; + for ( j = 0; j < numlegs; j++ ) { + edge = eg->nodes[i]->edges[j]; + vect = ABS(edge)-1; + if ( vect < info->numextern && info->legcouple[vect][numlegs] == 0 ) { + +// This cannot be. + + info->numtopo++; + return False; + } + } + } + } +#endif if ( ( info->flags & TOPOLOGIESONLY ) == 0 ) { info->numtopo++; return True; @@ -355,6 +513,8 @@ Bool ProcessTopology(EGraph *eg, void *ti) WORD *tail = tdia + tdia[1]; WORD *tend = term + *term; WORD *fill, *startfill; + Model *model = (Model *)info->currentModel; + MODEL *m = (MODEL *)info->currentMODEL; int i, j; int numlegs, vect, edge; @@ -374,6 +534,26 @@ Bool ProcessTopology(EGraph *eg, void *ti) *fill++ = 0; FILLFUN(fill) *fill++ = -SNUMBER; *fill++ = i+1; + if ( model != NULL && m != NULL ) { + if ( !eg->isExternal(i) ) { + WORD *afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill) + WORD *cfill = fill; *fill++ = 0; + int cpl = 2*eg->nodes[i]->extloop+eg->nodes[i]->deg-2; + *fill++ = SYMBOL; *fill++ = 4; + *fill++ = m->couplings[0]; + *fill++ = cpl; + *fill++ = 1; *fill++ = 1; *fill++ = 3; + *cfill = fill - cfill; + *afill = fill - afill; + if ( *afill == ARGHEAD+8 && afill[ARGHEAD+4] == 1 ) { + fill = afill; *fill++ = -SYMBOL; *fill++ = afill[ARGHEAD+3]; + } + } + else { + *fill++ = -SNUMBER; + *fill++ = 1; + } + } // // Now the momenta. // @@ -390,6 +570,116 @@ Bool ProcessTopology(EGraph *eg, void *ti) } startfill[1] = fill-startfill; } + if ( ( info->flags & WITHEDGES ) == WITHEDGES ) { + for ( i = 0; i < eg->nEdges; i++ ) { + int n1 = eg->edges[i]->nodes[0]; + int n2 = eg->edges[i]->nodes[1]; +// int l1 = eg->edges[i]->nlegs[0]; +// int l2 = eg->edges[i]->nlegs[1]; + startfill = fill; + *fill++ = EDGE; + *fill++ = 0; + FILLFUN(fill) + *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge +// + *fill++ = -VECTOR; + if ( i < info->numextern ) { // Look up in set of external momenta + *fill++ = SetElements[Sets[info->externalset].first+i]; + } + else { // Look up in set of internal momenta set + *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)]; + } +// + *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from + *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to + startfill[1] = fill - startfill; + } + } +// +// Block information +// + if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) { + for ( i = 0; i < eg->econn->nblocks; i++ ) { + startfill = fill; + *fill++ = BLOCK; + *fill++ = 0; + FILLFUN(fill); + *fill++ = -SNUMBER; + *fill++ = i+1; + *fill++ = -SNUMBER; + *fill++ = eg->econn->blocks[i].loop; +// +// Now we have to make a list of all nodes inside this block +// + int bnodes[GRCC_MAXNODES], k; + WORD *argfill = fill, *funfill; + *fill++ = 0; *fill++ = 0; FILLARG(fill) + *fill++ = 0; + for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0; + for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) { + bnodes[eg->econn->blocks[i].edges[k][0]] = 1; + bnodes[eg->econn->blocks[i].edges[k][1]] = 1; + } + for ( k = 0; k < GRCC_MAXNODES; k++ ) { + if ( bnodes[k] == 0 ) continue; +// +// Now we put the node inside this argument. +// + funfill = fill; + *fill++ = NODEFUNCTION; + *fill++ = 0; + FILLFUN(fill) + *fill++ = -SNUMBER; *fill++ = k+1; + numlegs = eg->nodes[k]->deg; + for ( j = 0; j < numlegs; j++ ) { + edge = eg->nodes[k]->edges[j]; + vect = ABS(edge)-1; + *fill++ = -VECTOR; + if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta + *fill++ = SetElements[Sets[info->externalset].first+vect]; + } + else { // Look up in set of internal momenta set + *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)]; + } + } + funfill[1] = fill-funfill; + } + *fill++ = 1; *fill++ = 1; *fill++ = 3; + *argfill = fill - argfill; + argfill[ARGHEAD] = argfill[0] - ARGHEAD; + startfill[1] = fill-startfill; + } + +// if ( eg->econn->narticuls > 0 ) { +// startfill = fill; +// *fill++ = BLOCK; +// *fill++ = 0; +// FILLFUN(fill); +// for ( i = 0; i < eg->econn->snodes; i++ ) { +// if ( eg->econn->articuls[i] != 0 ) { +// *fill++ = -SNUMBER; +// *fill++ = i+1; +// } +// } +// startfill[1] = fill-startfill; +// } + } + if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { + for ( i = 0; i < eg->econn->nopic; i++ ) { + startfill = fill; + *fill++ = ONEPI; + *fill++ = 0; + FILLFUN(fill); + *fill++ = -SNUMBER; + *fill++ = i+1; + *fill++ = -ONEPI; + for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) { + *fill++ = -SNUMBER; + *fill++ = eg->econn->opics[i].nodes[j]+1; + } + startfill[1] = fill-startfill; + } + } // // Topology counter. We have exaggerated a bit with the eye on the far future. // @@ -428,13 +718,15 @@ Bool ProcessTopology(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; +//MesPrint("<> %a",*newterm,newterm); + Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; info->numtopo++; - return True; + return False; } -// #] ProcessTopology : +// #] ProcessTopology : // #[ GenDiagrams : int GenDiagrams(PHEAD WORD *term, WORD level) @@ -447,7 +739,7 @@ int GenDiagrams(PHEAD WORD *term, WORD level) int babble = 0; // Later we may set this at the FORM code level TERMINFO info; WORD inset,outset,*coupl,setnum,optionnumber = 0; - int i, cpl[GRCC_MAXNCPLG]; + int i, j, cpl[GRCC_MAXNCPLG]; int ninitl, initlPart[GRCC_MAXLEGS], nfinal, finalPart[GRCC_MAXLEGS]; for ( i = 0; i < GRCC_MAXNCPLG; i++ ) cpl[i] = 0; // @@ -459,6 +751,7 @@ int GenDiagrams(PHEAD WORD *term, WORD level) info.diaoffset = AR.funoffset; info.externalset = term[info.diaoffset+FUNHEAD+7]; info.internalset = term[info.diaoffset+FUNHEAD+9]; + info.flags = 0; inset = term[info.diaoffset+FUNHEAD+3]; outset = term[info.diaoffset+FUNHEAD+5]; coupl = term + info.diaoffset + FUNHEAD + 10; @@ -471,7 +764,6 @@ int GenDiagrams(PHEAD WORD *term, WORD level) if ( term[info.diaoffset+1] > *coupl+FUNHEAD+10 ) optionnumber = term[info.diaoffset+*coupl+FUNHEAD+11]; } - setnum = term[info.diaoffset+FUNHEAD+1]; m = AC.models[SetElements[Sets[setnum].first]]; @@ -486,33 +778,59 @@ int GenDiagrams(PHEAD WORD *term, WORD level) opt = new Options(); - opt->setOutAG(ProcessDiagram, (void*) &info); - opt->setOutMG(ProcessTopology, (void*) &info); + opt->setOutAG(ProcessDiagram, &info); + opt->setOutMG(ProcessTopology, &info); +// opt->setEndMG(fendMG, &info); opt->values[GRCC_OPT_1PI] = ( optionnumber & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; opt->values[GRCC_OPT_NoTadpole] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; +// +// Next are snails: +// opt->values[GRCC_OPT_No1PtBlock] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; - opt->values[GRCC_OPT_No2PtL1PI] = ( optionnumber & WITHINSERTIONS ) == WITHINSERTIONS; +// + if ( ( optionnumber & WITHINSERTIONS ) == WITHINSERTIONS ) { + opt->values[GRCC_OPT_No2PtL1PI] = True; + opt->values[GRCC_OPT_NoAdj2PtV] = True; + opt->values[GRCC_OPT_No2PtL1PI] = True; + } + else { + opt->values[GRCC_OPT_NoAdj2PtV] = True; + } opt->values[GRCC_OPT_SymmInitial] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; - if ( ( optionnumber & TOPOLOGIESONLY ) == TOPOLOGIESONLY ) - info.flags |= TOPOLOGIESONLY; + opt->values[GRCC_OPT_SymmFinal] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; + +// opt->values[GRCC_OPT_Block] = ( optionnumber & WITHBLOCKS ) == WITHBLOCKS; opt->setOutputF(False,""); + opt->setOutputP(False,""); opt->printLevel(babble); +// opt->values[GRCC_OPT_Step] = GRCC_AGraph; + // Load the various arrays. ninitl = Sets[inset].last - Sets[inset].first; for ( i = 0; i < ninitl; i++ ) { x = SetElements[Sets[inset].first+i]; initlPart[i] = ConvertParticle(model,x); + info.legcouple[i] = m->vertices[numParticle(m,x)]->couplings; } nfinal = Sets[outset].last - Sets[outset].first; for ( i = 0; i < nfinal; i++ ) { x = SetElements[Sets[outset].first+i]; finalPart[i] = ConvertParticle(model,x); + info.legcouple[i+ninitl] = m->vertices[numParticle(m,x)]->couplings; } info.numextern = ninitl + nfinal; + for ( i = 2; i <= MAXLEGS; i++ ) { + if ( m->legcouple[i] == 1 ) { + for ( j = 0; j < info.numextern; j++ ) { + if ( info.legcouple[j][i] == 0 ) { info.flags |= CHECKEXTERN; goto Go_on; } + } + } + } +Go_on:; // // Now we have to sort out the coupling constants. // The argument at coupl can be of type -SNUMBER, -SYMBOL or generic @@ -529,11 +847,20 @@ int GenDiagrams(PHEAD WORD *term, WORD level) int nc = coupl[1]*2 + ninitl + nfinal - 2; int *scratch = (int *)Malloc1(nc*sizeof(int),"DistrN"); scratch[0] = -2; // indicating startup cq first call. - while ( DistrN(nc,cpl,m->ncouplings,scratch) ) { + + if ( ( info.flags & TOPOLOGIESONLY ) == 0 ) { + while ( DistrN(nc,cpl,m->ncouplings,scratch) ) { + proc = new Process(pid, model, opt, + ninitl, initlPart, nfinal, finalPart, cpl); + delete proc; + info.numtopo = 1; + } + } + else { + cpl[0] = nc; proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); + ninitl, initlPart, nfinal, finalPart, cpl); delete proc; - info.numtopo = 1; } M_free(scratch,"DistrN"); opt->end(); @@ -562,6 +889,9 @@ int GenDiagrams(PHEAD WORD *term, WORD level) t += 2; } } +/* + And now the generation: +*/ proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl); opt->end(); @@ -596,15 +926,11 @@ int processVertex(TOPOTYPE *TopoInf, int pointsremaining, int level) TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level]; TopoInf->clnum[TopoInf->ncl] = j; TopoInf->clext[TopoInf->ncl] = 0; - TopoInf->cmind[TopoInf->ncl] = 0; - TopoInf->cmaxd[TopoInf->ncl] = 0; TopoInf->ncl++; - MGraph *mgraph = new MGraph(1, - TopoInf->ncl, TopoInf->cldeg, - TopoInf->clnum, TopoInf->clext, - TopoInf->cmind, TopoInf->cmaxd, - TopoInf->opt); + MGraph *mgraph = new MGraph(1, TopoInf->ncl, TopoInf->cldeg, + TopoInf->clnum, TopoInf->clext, + TopoInf->cmind, TopoInf->cmaxd, TopoInf->opt); mgraph->generate(); @@ -663,6 +989,8 @@ int GenTopologies(PHEAD WORD *term, WORD level) 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]); @@ -674,19 +1002,27 @@ int GenTopologies(PHEAD WORD *term, WORD level) if ( t1 < tstop && t1[0] == -SNUMBER ) { if ( ( t1[1] & NONODES ) == NONODES ) info.flags |= NONODES; if ( ( t1[1] & WITHEDGES ) == WITHEDGES ) info.flags |= WITHEDGES; + if ( ( t1[1] & WITHBLOCKS ) == WITHBLOCKS ) info.flags |= WITHBLOCKS; + if ( ( t1[1] & WITHONEPI ) == WITHONEPI ) info.flags |= WITHONEPI; opt->values[GRCC_OPT_1PI] = ( t1[1] & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; - opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; +// opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; + opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOSNAILS ) == NOSNAILS; opt->values[GRCC_OPT_No1PtBlock] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; - opt->values[GRCC_OPT_No2PtL1PI] = ( t1[1] & WITHINSERTIONS ) == WITHINSERTIONS; + 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] & WITHSYMMETRIZE ) == WITHSYMMETRIZE; } info.numdia = 0; info.numtopo = 1; - opt->setOutAG(ProcessDiagram, (void*) &info); - opt->setOutMG(ProcessTopology, (void*) &info); - + 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 @@ -701,7 +1037,6 @@ int GenTopologies(PHEAD WORD *term, WORD level) } for ( i = 0; i < nlegs; i++ ) { TopoInf.cldeg[i] = 1; TopoInf.clnum[i] = 1; TopoInf.clext[i] = -1; - TopoInf.cmind[i] = 0; TopoInf.cmaxd[i] = 0; } int points = 2*nloops-2+nlegs; From ab36b49f50fd5c783fcc5a890305da0420ea4c84 Mon Sep 17 00:00:00 2001 From: jodavies Date: Tue, 18 Nov 2025 21:57:41 +0000 Subject: [PATCH 2/3] feat(grcc): implment qgraf-compatible keywords in FORM These pass a large number of comparisons with qgraf, using @tueda's "qgraf.frm" comparison procedures. --- sources/diawrap.cc | 142 +++++++++++++++++++++++++------------------- sources/ftypes.h | 44 ++++++++++---- sources/grcc.cc | 10 +++- sources/grccparam.h | 3 +- sources/startup.c | 43 ++++++++++---- 5 files changed, 154 insertions(+), 88 deletions(-) diff --git a/sources/diawrap.cc b/sources/diawrap.cc index c6d4ddb4..e9b0800f 100644 --- a/sources/diawrap.cc +++ b/sources/diawrap.cc @@ -245,7 +245,7 @@ void ProcessDiagram(EGraph *eg, void *ti) // // Now get the nodes // - if ( ( info->flags & NONODES ) == 0 ) { + if ( ( info->flags & WITHOUTNODES ) == 0 ) { for ( i = 0; i < eg->nNodes; i++ ) { // // node_(number,coupling,particle_1(momentum_1),...,particle_n(momentum_n)) @@ -395,7 +395,7 @@ void ProcessDiagram(EGraph *eg, void *ti) startfill[1] = fill-startfill; } } - if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { + if ( ( info->flags & WITHONEPISETS ) == WITHONEPISETS ) { for ( i = 0; i < eg->econn->nopic; i++ ) { startfill = fill; *fill++ = ONEPI; @@ -452,8 +452,6 @@ void ProcessDiagram(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; -// MesPrint("<> %a",newterm[0],newterm); - Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; } @@ -479,6 +477,9 @@ Bool ProcessTopology(EGraph *eg, void *ti) // return False; skip diagram generation (when asked for). // TERMINFO *info = (TERMINFO *)ti; + +// This seems to work properly. It was disabled before. +#define WITHEARLYVETO #ifdef WITHEARLYVETO if ( ( ( info->flags & CHECKEXTERN ) == CHECKEXTERN ) && info->currentMODEL != NULL ) { int i, j; @@ -500,6 +501,7 @@ Bool ProcessTopology(EGraph *eg, void *ti) } } #endif + if ( ( info->flags & TOPOLOGIESONLY ) == 0 ) { info->numtopo++; return True; @@ -664,7 +666,7 @@ Bool ProcessTopology(EGraph *eg, void *ti) // startfill[1] = fill-startfill; // } } - if ( ( info->flags & WITHONEPI ) == WITHONEPI ) { + if ( ( info->flags & WITHONEPISETS ) == WITHONEPISETS ) { for ( i = 0; i < eg->econn->nopic; i++ ) { startfill = fill; *fill++ = ONEPI; @@ -718,8 +720,6 @@ Bool ProcessTopology(EGraph *eg, void *ti) *newterm = fill - newterm; AT.WorkPointer = fill; -//MesPrint("<> %a",*newterm,newterm); - Generator(BHEAD newterm,info->level); AT.WorkPointer = oldworkpointer; info->numtopo++; @@ -727,6 +727,32 @@ Bool ProcessTopology(EGraph *eg, void *ti) } // #] ProcessTopology : +// #[ SetDualOpts : +void SetDualOpts(int *opt, const WORD num, const int key, const char* key_name, + const int dual, const char* dual_name, const int val, const int dval) { + + if ( ( num & key ) == key ) { + if ( ( num & dual ) == dual ) { + MLOCK(ErrorMessageLock); + MesPrint("&Conflicting diagram filters: %s and %s.", key_name, dual_name); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } + else { + *opt = val; + } + } + else { + if ( ( num & dual ) == dual ) { + *opt = dval; + } + else { + // The default value is always 0. + *opt = 0; + } + } +} +// #] SetDualOpts : // #[ GenDiagrams : int GenDiagrams(PHEAD WORD *term, WORD level) @@ -780,48 +806,46 @@ int GenDiagrams(PHEAD WORD *term, WORD level) opt->setOutAG(ProcessDiagram, &info); opt->setOutMG(ProcessTopology, &info); -// opt->setEndMG(fendMG, &info); - opt->values[GRCC_OPT_1PI] = ( optionnumber & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; - opt->values[GRCC_OPT_NoTadpole] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; -// -// Next are snails: -// - opt->values[GRCC_OPT_No1PtBlock] = ( optionnumber & NOTADPOLES ) == NOTADPOLES; -// - if ( ( optionnumber & WITHINSERTIONS ) == WITHINSERTIONS ) { - opt->values[GRCC_OPT_No2PtL1PI] = True; - opt->values[GRCC_OPT_NoAdj2PtV] = True; - opt->values[GRCC_OPT_No2PtL1PI] = True; - } - else { - opt->values[GRCC_OPT_NoAdj2PtV] = True; - } - opt->values[GRCC_OPT_SymmInitial] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; - opt->values[GRCC_OPT_SymmFinal] = ( optionnumber & WITHSYMMETRIZE ) == WITHSYMMETRIZE; - -// opt->values[GRCC_OPT_Block] = ( optionnumber & WITHBLOCKS ) == WITHBLOCKS; + opt->values[GRCC_OPT_SymmInitial] = ( optionnumber & WITHSYMMETRIZEI ) == WITHSYMMETRIZEI; + opt->values[GRCC_OPT_SymmFinal] = ( optionnumber & WITHSYMMETRIZEF ) == WITHSYMMETRIZEF; + + // WITHBLOCKS controls output formatting. We could introduce an extra filtering option + // corresponding to GRCC_OPT_Block, which is somewhat like Qgraf "onevi" but not quite + // the same currently. + //opt->values[GRCC_OPT_Block] = ; + + // Now the "qgraf-compatible filtering options": + int qgopt[GRCC_QGRAF_OPT_Size]; + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_ONEPI], optionnumber,ONEPARTI, "ONEPI_", ONEPARTR, "ONEPR_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_ONSHELL], optionnumber,ONSHELL, "ONSHELL_", OFFSHELL, "OFFSHELL_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOSIGMA], optionnumber,NOSIGMA, "NOSIGMA_", SIGMA, "SIGMA_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOSNAIL], optionnumber,NOSNAIL, "NOSNAIL_", SNAIL, "SNAIL_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOTADPOLE],optionnumber,NOTADPOLE,"NOTADPOLE_",TADPOLE , "TADPOLE_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_SIMPLE], optionnumber,SIMPLE, "SIMPLE_", NOTSIMPLE,"NOTSIMPLE_",1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_BIPART], optionnumber,BIPART, "BIPART_", NONBIPART,"NONBIPART_",1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_CYCLI], optionnumber,CYCLI, "CYCLI_", CYCLR, "CYCLR_", 1,-1); + SetDualOpts(&qgopt[GRCC_QGRAF_OPT_FLOOP], optionnumber,FLOOP, "FLOOP_", NOTFLOOP, "NOTFLOOP_", 1,-1); + // Now set the options internally: + opt->setQGrafOpt(qgopt); opt->setOutputF(False,""); opt->setOutputP(False,""); opt->printLevel(babble); -// opt->values[GRCC_OPT_Step] = GRCC_AGraph; - // Load the various arrays. - - ninitl = Sets[inset].last - Sets[inset].first; - for ( i = 0; i < ninitl; i++ ) { - x = SetElements[Sets[inset].first+i]; - initlPart[i] = ConvertParticle(model,x); + ninitl = Sets[inset].last - Sets[inset].first; + for ( i = 0; i < ninitl; i++ ) { + x = SetElements[Sets[inset].first+i]; + initlPart[i] = ConvertParticle(model,x); info.legcouple[i] = m->vertices[numParticle(m,x)]->couplings; - } - nfinal = Sets[outset].last - Sets[outset].first; - for ( i = 0; i < nfinal; i++ ) { - x = SetElements[Sets[outset].first+i]; - finalPart[i] = ConvertParticle(model,x); + } + nfinal = Sets[outset].last - Sets[outset].first; + for ( i = 0; i < nfinal; i++ ) { + x = SetElements[Sets[outset].first+i]; + finalPart[i] = ConvertParticle(model,x); info.legcouple[i+ninitl] = m->vertices[numParticle(m,x)]->couplings; - } + } info.numextern = ninitl + nfinal; for ( i = 2; i <= MAXLEGS; i++ ) { if ( m->legcouple[i] == 1 ) { @@ -850,16 +874,14 @@ Go_on:; if ( ( info.flags & TOPOLOGIESONLY ) == 0 ) { while ( DistrN(nc,cpl,m->ncouplings,scratch) ) { - proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); + proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl); delete proc; info.numtopo = 1; } } else { cpl[0] = nc; - proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); + proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl); delete proc; } M_free(scratch,"DistrN"); @@ -892,8 +914,7 @@ Go_on:; /* And now the generation: */ - proc = new Process(pid, model, opt, - ninitl, initlPart, nfinal, finalPart, cpl); + proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl); opt->end(); delete proc; delete opt; @@ -1000,22 +1021,23 @@ int GenTopologies(PHEAD WORD *term, WORD level) info.flags |= TOPOLOGIESONLY; // this is the topologies_ function after all. if ( t1 < tstop && t1[0] == -SNUMBER ) { - if ( ( t1[1] & NONODES ) == NONODES ) info.flags |= NONODES; + 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] & WITHONEPI ) == WITHONEPI ) info.flags |= WITHONEPI; - opt->values[GRCC_OPT_1PI] = ( t1[1] & ONEPARTICLEIRREDUCIBLE ) == ONEPARTICLEIRREDUCIBLE; -// opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; - opt->values[GRCC_OPT_NoTadpole] = ( t1[1] & NOSNAILS ) == NOSNAILS; - opt->values[GRCC_OPT_No1PtBlock] = ( t1[1] & NOTADPOLES ) == NOTADPOLES; - 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] & WITHSYMMETRIZE ) == WITHSYMMETRIZE; + 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; diff --git a/sources/ftypes.h b/sources/ftypes.h index 80246c44..be1a8439 100644 --- a/sources/ftypes.h +++ b/sources/ftypes.h @@ -1102,16 +1102,36 @@ typedef int (*TFUN1)(); #define DENSETABLE 1 #define SPARSETABLE 0 -#define ONEPARTICLEIRREDUCIBLE 1 -#define WITHINSERTIONS 2 -#define NOTADPOLES 4 -#define WITHSYMMETRIZE 8 -#define TOPOLOGIESONLY 16 -#define NONODES 32 -#define WITHEDGES 64 -#define CHECKEXTERN 128 -#define WITHBLOCKS 256 -#define WITHONEPI 512 -#define NOSNAILS 1024 -#define NOEXTSELF 2048 +// Diagram generator flags. They should be powers of two, since they are added +// to pass to diagrams and masked in diawrap.cc to check the flags. +// We use a "stringified" version of these when defining the FORM preprocessor +// variables, so can't neatly use "1<<10" etc here. +#define TOPOLOGIESONLY 1 +#define WITHOUTNODES 2 +#define WITHEDGES 4 +#define WITHBLOCKS 8 +#define WITHONEPISETS 16 +#define WITHSYMMETRIZEI 32 +#define WITHSYMMETRIZEF 64 +// This is not an "option" preprocessor var but is set by "external" particle definitions +#define CHECKEXTERN 128 +// The "qgraf compatible" filtering flags: +#define ONEPARTI 256 +#define ONEPARTR 512 +#define ONSHELL 1024 +#define OFFSHELL 2048 +#define NOSIGMA 4096 +#define SIGMA 8192 +#define NOSNAIL 16384 +#define SNAIL 32768 +#define NOTADPOLE 65536 +#define TADPOLE 131072 +#define SIMPLE 262144 +#define NOTSIMPLE 524288 +#define BIPART 1048576 +#define NONBIPART 2097152 +#define CYCLI 4194304 +#define CYCLR 8388608 +#define FLOOP 16777216 +#define NOTFLOOP 33554432 diff --git a/sources/grcc.cc b/sources/grcc.cc index 18089d89..8272f86e 100644 --- a/sources/grcc.cc +++ b/sources/grcc.cc @@ -7242,16 +7242,22 @@ Bool EGraph::optQGrafA(Options *opt) printf("optQGrafA: %8ld\n", mId); econn->print(); #endif + Bool retval = True; if (opt->qgopt[GRCC_QGRAF_OPT_FLOOP] != 0) { for (int fl=0; fl < nFlines; fl++) { if (flines[fl]->ftype == FL_Closed) { if (flines[fl]->nlist % 2 != 0) { - return False; + retval = False; + break; } } } + // `notfloop_' is the dual of `floop_', so flip the decision: + if (opt->qgopt[GRCC_QGRAF_OPT_FLOOP] == -1) { + retval = (retval == True ? False : True); + } } - return True; + return retval; } //-------------------------------------------------------------- diff --git a/sources/grccparam.h b/sources/grccparam.h index e769b0b5..88d03e5d 100644 --- a/sources/grccparam.h +++ b/sources/grccparam.h @@ -129,7 +129,8 @@ typedef struct { #define GRCC_QGRAF_OPT_BIPART 6 #define GRCC_QGRAF_OPT_CYCLI 7 #define GRCC_QGRAF_OPT_FLOOP 8 -#define GRCC_QGRAF_OPT_TOPOL 9 +//TODO what does this do? Does it work? +//#define GRCC_QGRAF_OPT_TOPOL 9 #ifdef GRCC_QGRAF_OPT_TOPOL #define GRCC_QGRAF_OPT_Size 10 diff --git a/sources/startup.c b/sources/startup.c index 7de4474b..1b39678f 100644 --- a/sources/startup.c +++ b/sources/startup.c @@ -1205,20 +1205,37 @@ void StartVariables(void) PutPreVar((UBYTE *)"keepright_",(UBYTE *)("0"),0,0); PutPreVar((UBYTE *)"SYSTEMERROR_",(UBYTE *)("0"),0,0); /* - Next are a few 'constants' for diagram generation + Next are the flags to control diagram generation filters */ - PutPreVar((UBYTE *)"ONEPI_",(UBYTE *)("1"),0,0); - PutPreVar((UBYTE *)"WITHOUTINSERTIONS_",(UBYTE *)("2"),0,0); - PutPreVar((UBYTE *)"NOTADPOLES_",(UBYTE *)("4"),0,0); - PutPreVar((UBYTE *)"SYMMETRIZE_",(UBYTE *)("8"),0,0); - PutPreVar((UBYTE *)"TOPOLOGIESONLY_",(UBYTE *)("16"),0,0); - PutPreVar((UBYTE *)"NONODES_",(UBYTE *)("32"),0,0); - PutPreVar((UBYTE *)"WITHEDGES_",(UBYTE *)("64"),0,0); -/* Note that CHECKEXTERN is 128 */ - PutPreVar((UBYTE *)"WITHBLOCKS_",(UBYTE *)("256"),0,0); - PutPreVar((UBYTE *)"WITHONEPISETS_",(UBYTE *)("512"),0,0); - PutPreVar((UBYTE *)"NOSNAILS_",(UBYTE *)("1024"),0,0); - PutPreVar((UBYTE *)"NOEXTSELF_",(UBYTE *)("2048"),0,0); + #define STR2(x) #x + #define STR(x) STR2(x) + PutPreVar((UBYTE *)"TOPOLOGIESONLY_" ,(UBYTE*)(STR(TOPOLOGIESONLY)) ,0,0); + PutPreVar((UBYTE *)"WITHOUTNODES_" ,(UBYTE*)(STR(WITHOUTNODES)) ,0,0); + PutPreVar((UBYTE *)"WITHEDGES_" ,(UBYTE*)(STR(WITHEDGES)) ,0,0); + PutPreVar((UBYTE *)"WITHBLOCKS_" ,(UBYTE*)(STR(WITHBLOCKS)) ,0,0); + PutPreVar((UBYTE *)"WITHONEPISETS_" ,(UBYTE*)(STR(WITHONEPISETS)) ,0,0); + PutPreVar((UBYTE *)"WITHSYMMETRIZEI_",(UBYTE*)(STR(WITHSYMMETRIZEI)),0,0); + PutPreVar((UBYTE *)"WITHSYMMETRIZEF_",(UBYTE*)(STR(WITHSYMMETRIZEF)),0,0); + // This is not an "option" preprocessor var but is set by "external" particle definitions +// PutPreVar((UBYTE *)"CHECKEXTERN_" ,(UBYTE*)(STR(CHECKEXTERN)) ,0,0); + PutPreVar((UBYTE *)"ONEPI_" ,(UBYTE*)(STR(ONEPARTI)) ,0,0); + PutPreVar((UBYTE *)"ONEPR_" ,(UBYTE*)(STR(ONEPARTR)) ,0,0); + PutPreVar((UBYTE *)"ONSHELL_" ,(UBYTE*)(STR(ONSHELL)) ,0,0); + PutPreVar((UBYTE *)"OFFSHELL_" ,(UBYTE*)(STR(OFFSHELL)) ,0,0); + PutPreVar((UBYTE *)"NOSIGMA_" ,(UBYTE*)(STR(NOSIGMA)) ,0,0); + PutPreVar((UBYTE *)"SIGMA_" ,(UBYTE*)(STR(SIGMA)) ,0,0); + PutPreVar((UBYTE *)"NOSNAIL_" ,(UBYTE*)(STR(NOSNAIL)) ,0,0); + PutPreVar((UBYTE *)"SNAIL_" ,(UBYTE*)(STR(SNAIL)) ,0,0); + PutPreVar((UBYTE *)"NOTADPOLE_" ,(UBYTE*)(STR(NOTADPOLE)) ,0,0); + PutPreVar((UBYTE *)"TADPOLE_" ,(UBYTE*)(STR(TADPOLE)) ,0,0); + PutPreVar((UBYTE *)"SIMPLE_" ,(UBYTE*)(STR(SIMPLE)) ,0,0); + PutPreVar((UBYTE *)"NOTSIMPLE_" ,(UBYTE*)(STR(NOTSIMPLE)) ,0,0); + PutPreVar((UBYTE *)"BIPART_" ,(UBYTE*)(STR(BIPART)) ,0,0); + PutPreVar((UBYTE *)"NONBIPART_" ,(UBYTE*)(STR(NONBIPART)) ,0,0); + PutPreVar((UBYTE *)"CYCLI_" ,(UBYTE*)(STR(CYCLI)) ,0,0); + PutPreVar((UBYTE *)"CYCLR_" ,(UBYTE*)(STR(CYCLR)) ,0,0); + PutPreVar((UBYTE *)"FLOOP_" ,(UBYTE*)(STR(FLOOP)) ,0,0); + PutPreVar((UBYTE *)"NOTFLOOP_" ,(UBYTE*)(STR(NOTFLOOP)) ,0,0); { char buf[41]; /* up to 128-bit */ From 9dc826bb96590b1751654c13a64d2458e20da54d Mon Sep 17 00:00:00 2001 From: Takahiro Ueda Date: Mon, 24 Nov 2025 22:40:42 +0900 Subject: [PATCH 3/3] test: diagram generator comparison tests Add FORM routines for comparing the diagram generator with QGRAF and FeynGraph: http://cefema-gt.tecnico.ulisboa.pt/~paulo/qgraf.html https://jens-braun.github.io/FeynGraph/ --- check/diagrams/compare.frm | 1060 ++++++++++++++++++++++++++++++++++++ 1 file changed, 1060 insertions(+) create mode 100644 check/diagrams/compare.frm diff --git a/check/diagrams/compare.frm b/check/diagrams/compare.frm new file mode 100644 index 00000000..766c43cf --- /dev/null +++ b/check/diagrams/compare.frm @@ -0,0 +1,1060 @@ +* Compares the diagram generator with others. +#ifndef `TEST' + #message Use -D TEST=XXX + #terminate +#else + #include `NAME_' # `TEST' +#endif +.end + +*--#[ diagram_compare_include : + +*#define Verbose +*#define CacheRead +*#define CacheWrite +#ifndef `CacheDir' + #define CacheDir "." +#endif + +#ifdef `TESTFILEDIR' + #redefine CacheRead + #redefine CacheDir "`TESTFILEDIR'/cache" +#endif + +#define MaxLegs "9" +#define MaxPropagators "19" + +#define CurrentModel "" +#define Bosons "" +#define AntiParticles "" + +#define DiagramGenerator "" +#define DiagramGeneratorVersion "" + +#define QgrafVersion "" +#define FeynGraphVersion "" + +#define QgrafStyleFile "qgraf.sty" +#define QgrafModelFile "qgraf.mdl" +#define QgrafInputFile "qgraf.dat" +#define QgrafOutputFile "qgraf.out" +#define FeynGraphTemplateFile "feyngraph.jinja" +#define FeynGraphInputFile "feyngraph.toml" +#define FeynGraphOutputFile "diagrams.out" + +Format nospaces; +Format 255; + +* Momentum convention: all incoming +* incoming particles: p1, p2, ... +* outgoing particles: -q1, -q2, ... +* internal particles: k1, k2, ... + +S n1,n2,x1,x2; +CF topo,node,counter,vx,vxs(symmetric),acc,replace; +V p1,...,p`MaxLegs',q1,...,q`MaxLegs',k1,...,k`MaxPropagators',l1,...,l`MaxPropagators'; +Set kk:k1,...,k`MaxPropagators'; +Set ll:l1,...,l`MaxPropagators'; + +#OpenDictionary wildmom + #do i = 1, `MaxPropagators' + #add k`i': "k`i'?$k`i'" + #enddo +#CloseDictionary + +#do i = 1, `MaxPropagators' + #$k`i' = 0; +#enddo + +Model TEMP; + Particle f1; + Particle f2; +EndModel; + +*--#[ SetupPhi3Model : + +** +* Prepares the PHI3 model. +* +#procedure SetupPhi3Model() + #define filename "`QgrafModelFile'" + #create <`filename'> + #write <`filename'> "[phi,phi,+]" + #write <`filename'> "[phi,phi,phi]" + #close <`filename'> + + Model PHI3; + Particle phi,1; + Vertex phi,phi,phi:g; + EndModel; + + #redefine CurrentModel "PHI3" + #redefine Bosons "phi" + #redefine AntiParticles "" +#endprocedure + +*--#] SetupPhi3Model : +*--#[ SetupPhi4Model : + +** +* Prepares the PHI4 model. +* +#procedure SetupPhi4Model() + #define filename "`QgrafModelFile'" + #create <`filename'> + #write <`filename'> "[phi,phi,+]" + #write <`filename'> "[phi,phi,phi,phi]" + #close <`filename'> + + Model PHI4; + Particle phi,1; + Vertex phi,phi,phi,phi:g^2; + EndModel; + + #redefine CurrentModel "PHI4" + #redefine Bosons "phi" + #redefine AntiParticles "" +#endprocedure + +*--#] SetupPhi4Model : +*--#[ SetupQcdModel : + +** +* Prepares the QCD model. +* +#procedure SetupQcdModel() + #define filename "`QgrafModelFile'" + #create <`filename'> + #write <`filename'> "[qua,QUA,-]" + #write <`filename'> "[glu,glu,+]" + #write <`filename'> "[gho,GHO,-]" + #write <`filename'> "[QUA,qua,glu]" + #write <`filename'> "[glu,glu,glu]" + #write <`filename'> "[glu,glu,glu,glu]" + #write <`filename'> "[GHO,gho,glu]" + #close <`filename'> + + Model QCD; + Particle qua,QUA,-2; + Particle gho,GHO,-1; + Particle glu,+3; + Vertex QUA,qua,glu:g; + Vertex GHO,gho,glu:g; + Vertex glu,glu,glu:g; + Vertex glu,glu,glu,glu:g^2; + EndModel; + + #redefine CurrentModel "QCD" + #redefine Bosons "glu" + #redefine AntiParticles "QUA,GHO" +#endprocedure + +*--#] SetupQcdModel : +*--#[ MakeDiagrams : + +** +* Manipulates diagrams for further manipulation. +* +* Parameters +* ---------- +* numIn: positive integer +* The number of incoming particles. +* numOut: positive integer +* The number of outgoing particles. +* +* Example Input +* ------------- +* topo_(1) +* *node_(1,1,QUA(-p1)) +* *node_(2,1,qua(-q1)) +* *node_(3,g,QUA(k2),qua(p1),glu(k1)) +* *node_(4,g,QUA(q1),qua(-k2),glu(-k1)) +* +* vx(QUA(-2),qua(1),glu(2)) +* *vx(QUA(1),qua(-1),glu(2)) +* +* Example Output +* -------------- +* node(QUA(q1),qua(-k2),glu(-k1)) +* *node(QUA(k2),qua(p1),glu(k1)) +* +#procedure MakeDiagrams(numIn,numOut) +* for diagrams_ + multiply replace_(topo_,dummy_); + id node_(n1?,n2?,?a) = node_(?a); + id node_(f1?(p1?)) = 1; + multiply replace_(node_,node); + +* for QGRAF +* Incoming fields have indices -1, -3, -5, ..., - 2j+1, ... +* Outgoing fields have indices -2, -4, -6, ..., - 2j, ... +* See the manual of QGRAF section 5.1 "The indices". + #do i = 1, `numIn' + repeat id vx(?a,f1?(1-2*`i'),?b) = vx(?a,f1(p`i'),?b); + #enddo + #do i = 1, `numOut' + repeat id vx(?a,f1?(-2*`i'),?b) = vx(?a,f1(q`i'),?b); + #enddo + repeat id vx(?a,f1?(-1),?b) = vx(?a,f1(p1),?b); + repeat id vx(?a,f1?(-2),?b) = vx(?a,f1(q1),?b); + multiply counter(1); + repeat id vx(?a,f1?(n1?),?b) * vx(?c,f2?(n1?),?d) * counter(n2?) + = vx(?a,f1(kk[n2]),?b) * vx(?c,f2(-kk[n2]),?d) * counter(n2+1); + repeat id vx(?a,f1?(n1?),?b,f2?(n1?),?c) * counter(n2?) + = vx(?a,f1(kk[n2]),?b,f2(-kk[n2]),?c) * counter(n2+1); + multiply replace_(counter,dummy_); + multiply replace_(vx,node); + +#ifdef `Verbose' + P +sss; +#endif + .sort:mkdiag; +#endprocedure + +*--#] MakeDiagrams : +*--#[ CanonicalizeDiagrams : + +** +* Canonicalizes diagrams. +* +* This procedure attaches a unique topology ID and reassigns +* the loop momenta k1, k2, ... in a canonical way. +* +* Example Input +* ------------- +* node(QUA(q1),qua(-k2),glu(-k1)) +* *node(QUA(k2),qua(p1),glu(k1)) +* +* Example Output +* -------------- +* topo(1) +* *node(QUA(q1),qua(-k1),glu(k2)) +* *node(QUA(k1),qua(p1),glu(-k2)) +* +#procedure CanonicalizeDiagrams() + #define numactiveexprs "`numactiveexprs_'" + #define activeexprnames "`activeexprnames_'" + +* Label each term uniquely. + #$counter = 0; + #define i "0" + #do e = {`activeexprnames',} + #ifdef `e' + #redefine i "{`i'+1}" + #$counter`i' = 0; + inexpression `e'; + $counter`i' = $counter`i' + 1; + $counter = $counter`i'; + multiply counter($counter`i'); + goto end; + endinexpression; + #endif + #enddo + #undefine i +label end; + B node; + ModuleOption maximum, $counter; + #do i = 1, `numactiveexprs' + ModuleOption local, $counter`i'; + #enddo + Moduleoption inparallel; + .sort:canondiag-1; + +* Extract vertex structure into acc(vxs(v1,...)*...). + Keep Brackets; + id node(?a) = node(?a) * vxs(?a); + argument vxs; + id f1?(p1?) = p1; + endargument; + multiply acc(1); + multiply replace_(vxs,vx); + id vx(?a) = acc(vxs(?a)); + repeat id acc(x1?) * acc(x2?) = acc(x1 * x2); + B acc; + .sort:canondiag-2; + + #do e = {`activeexprnames',} + #ifdef `e' + CreateSpectator [`e'.spec], "`e'.`PID_'.spec"; + #endif + #enddo + + #$numtopo = 0; + + #do loop = 1, 1 +* Find a new topology. + Keep Brackets; + #$found = 0; + if ($found == 0); + if (match(acc(x1?$x1))); + $found = 1; + $numtopo = $numtopo + 1; + $topo = $x1; + inside $x1; + $numv = count_(vxs,1); + endinside; + endif; + endif; + B acc; + ModuleOption noparallel; + .sort:canondiag-3-{`$numtopo'+1}; + #if `$found' +* Try to match the topology with the current term in all possible ways. + Keep Brackets; + #$found = 0; + $replaced = 0; + argument acc; +#UseDictionary wildmom($) + if ((count(vxs,1) == `$numv') && match(`$topo')); +#CloseDictionary + $found = 1; + $replaced = 1; + $t = term_; + inside $t; + id all, $topo * replace_(,...,) + = $topo * *...*; + multiply replace_(vxs,dummy_); + id replace(-k1?vector_,k2?) = replace(k1,-k2); + chainin replace; + multiply replace_(,...,); + endinside; + multiply replace_(vxs,dummy_); + multiply topo(`$numtopo'); + multiply $t; + endif; + endargument; + if ($replaced); + id acc(x1?) = x1; + endif; + ModuleOption maximum, $found; + ModuleOption local, $t, $replaced, <$k1>,...,<$k`MaxPropagators'>; + .sort:canondiag-4-`$numtopo'; + #if `$found' +* Perform all the generated replacements. + #redefine loop "0" + id ifnomatch-> notreplaced, replace(?a) = replace_(?a); + #do e = {`activeexprnames',} + #ifdef `e' + inexpression `e'; + tospectator [`e'.spec]; + endinexpression; + #endif + #enddo +label notreplaced; + B acc; + .sort:canondiag-5-`$numtopo'; + #endif + #endif + #enddo + +* Restore expressions from the spectators. + #do e = {`activeexprnames',} + #ifdef `e' + CopySpectator [`e'.copy] = [`e'.spec]; + #exchange `e',[`e'.copy] + #endif + #enddo + .sort:canondiag-6; + #do e = {`activeexprnames',} + #ifdef `e' + Drop [`e'.copy]; + EmptySpectator [`e'.spec]; + #endif + #enddo + +* Drop duplicate terms (keep the first one). + #define maxunrolling "100" + #define unrolling + #define i "0" + #do k = 1, `$counter' + #if {`k'+`maxunrolling'-1} <= `$counter' + #redefine unrolling "`maxunrolling'" + #else + #redefine unrolling "{`$counter'-`k'+1}" + #endif + #do j = 0, {`unrolling'-1} + if (match(counter({`k'+`j'}))); + #redefine i "0" + #do e = {`activeexprnames',} + #ifdef `e' + #redefine i "{`i'+1}" + #$found`i'x`j' = 0; + inexpression `e'; + if ($found`i'x`j'); + discard; + else; + $found`i'x`j' = 1; + multiply replace_(counter,dummy_); + tospectator [`e'.spec]; + endif; + endinexpression; + #endif + #enddo + endif; + #enddo + #do i = 1, `numactiveexprs' + #do j = 0, {`unrolling'-1} + ModuleOption local, $found`i'x`j'; + #enddo + #enddo + ModuleOption inparallel; + .sort:canondiag-7-`k'; + #redefine k "{`k'+`unrolling'-1}" + #enddo + #undefine i + #if `$counter' == 0 + .sort:canondiag-7; + #endif + +* Restore expressions from the spectators. + #do e = {`activeexprnames',} + #ifdef `e' + CopySpectator [`e'.copy] = [`e'.spec]; + #exchange `e',[`e'.copy] + #endif + #enddo + .sort:canondiag-8; + #do e = {`activeexprnames',} + #ifdef `e' + Drop [`e'.copy]; + RemoveSpectator [`e'.spec]; + #endif + #enddo + +* Reorder identical boson fields at each vertex to align the momenta. + #do b = {`Bosons',} + #ifdef `b' +* Assume that they appear in adjacent positions. + repeat id node(?a,`b'(?b),`b'(?c),?d) = node(?a,`b'(?b,?c),?d); + multiply replace_(`b',vxs); + multiply replace_(vxs,`b'); + repeat id node(?a,`b'(p1?,p2?,?b),?c) = node(?a,`b'(p1),`b'(p2,?b),?c); + #endif + #enddo + +* Align the loop momenta with the flow of particles. + #ifdef `AntiParticles' + #do i = 1, `MaxPropagators' + if (match(node(?a,f1?{`AntiParticles'}(-k`i'),?b))); + multiply replace_(k`i',-k`i'); + endif; + #enddo + #endif + + .sort:canondiag-9; +#endprocedure + +*--#] CanonicalizeDiagrams : +*--#[ CompareDiagrams : + +** +* Compares diagrams in two expressions. +* +* The common diagrams are extracted from F1 and F2, and stored in F0. +* Diagrams must consist of "topo" and "node" functions. +* +* Parameters +* ---------- +* F1, F2: expression [in/out] +* The expressions to be compared. Common diagrams will be removed from +* these expressions. +* F0: expression [out] +* The extracted common diagrams. +* +#procedure CompareDiagrams(F1,F2,F0) + #call CanonicalizeDiagrams() + + B+ topo,node; + .sort:compdiag-1; + Keep Brackets; + L `F0' = 0 + #do t = `F1'; + #$t = `t'; + #inside $t; + dropcoefficient; + #endinside; + #$t1 = `F1'[`$t']; + #$t2 = `F2'[`$t']; + #$tt = $t1 - $t2; + #if termsin($tt) == 0 + + (`t') + #endif + #enddo + ; + .sort + L `F1' = `F1' - `F0'; + L `F2' = `F2' - `F0'; +#ifdef `Verbose' + P +sss; +#endif + .sort:compdiag-3; +#endprocedure + +*--#] CompareDiagrams : +*--#[ DoComparison : + +** +* Performs comparison. +* +* Example +* ------- +* #call DoComparison(qcd,in=qua,out=qua,loops=1,form_options=`NOSNAIL_',qgraf_options=nosnail) +* +#procedure DoComparison(model,?a) + #switch `DiagramGenerator' + #case qgraf + #define DiagramGeneratorName "QGRAF" + #break + #case feyngraph + #define DiagramGeneratorName "FeynGraph" + #break + #default + #ifndef `DiagramGenerator' + #message DiagramGenerator not defined + #else + #message Unknown DiagramGenerator: `DiagramGenerator' + #endif + #terminate + #endswitch + + #define in "" + #define out "" + #define loops "" + #define formoptions "" + #define qgrafoptions "" + #define feyngraphoptions "" + + #define next "" + #do a = {`?a',} + #ifdef `a' + #if "`keepleft_(`a',3)'" == "in=" + #redefine next "in" + #redefine a "`takeleft_(`a',3)'" + #elseif "`keepleft_(`a',4)'" == "out=" + #redefine next "out" + #redefine a "`takeleft_(`a',4)'" + #elseif "`keepleft_(`a',6)'" == "loops=" + #redefine next "loops" + #redefine a "`takeleft_(`a',6)'" + #elseif "`keepleft_(`a',13)'" == "form_options=" + #redefine next "formoptions" + #redefine a "`takeleft_(`a',13)'" + #elseif "`keepleft_(`a',14)'" == "qgraf_options=" + #redefine next "qgrafoptions" + #redefine a "`takeleft_(`a',14)'" + #elseif "`keepleft_(`a',18)'" == "feyngraph_options=" + #redefine next "feyngraphoptions" + #redefine a "`takeleft_(`a',18)'" + #endif + #ifdef `a' + #ifdef `next' + #ifdef ``next'' + #redefine `next' "``next'',`a'" + #else + #redefine `next' "`a'" + #endif + #endif + #endif + #endif + #enddo + + #ifdef `formoptions' + #redefine formoptions "{`formoptions'}" + #else + #redefine formoptions "0" + #endif + + #define numIn "0" + #define numOut "0" + + #do a = {`in',} + #ifdef `a' + #redefine numIn "{`numIn'+1}" + #endif + #enddo + + #do a = {`out',} + #ifdef `a' + #redefine numOut "{`numOut'+1}" + #endif + #enddo + + #define momenta "" + + #do i = 1, `numIn' + #redefine momenta "`momenta',p`i'" + #enddo + + #do i = 1, `numOut' + #redefine momenta "`momenta',q`i'" + #enddo + + #define cachekey "`DiagramGenerator'-`model'" + #do b = {in,out,loops,formoptions,`DiagramGenerator'options} + #redefine cachekey "`cachekey'-" + #ifdef ``b'' + #do a = {``b'',} + #ifdef `a' + #redefine cachekey "`cachekey'`a'" + #endif + #enddo + #else + #redefine cachekey "`cachekey'none" + #endif + #enddo + + #define cachefile "`CacheDir'/`cachekey'.dat" + #define cachetmpfile "`cachekey'.dat.`PID_'.tmp" + #define cachetmpfile2 "`CacheDir'/`cachekey'.dat2.`PID_'.tmp" + +* Set up the model. + + #call Setup`toupper_(`keepleft_(`model',1)')'`takeleft_(`model',1)'Model() + +* Generate diagrams. + + L F1 = diagrams_(`CurrentModel',{`in'},{`out'},{`momenta'},kk,`loops',`formoptions'); +#ifdef `Verbose' + P +sss; +#endif + .sort:gen-diag; + #$nF1 = termsin_(F1); + +* Check if the cache exists. + + #ifdef `CacheRead' + #$cacheexists = + #pipe if [ -f "`cachefile'" ]; then echo 1; else echo 0; fi + ; + #else + #$cacheexists = 0; + #endif + + #if `$cacheexists' == 1 +* If the cache exists, then compare the generated diagrams with the cache, +* which has already been confirmed to be consistent with the other generator. + + #message Read cache file: `cachefile' + + Skip; + + L F2 = + #include `cachefile' + ; + .sort:load-cache; + Drop; + L oldF1 = F1; + L oldF2 = F2; + .sort; + Drop; + L F0 = 0; + L F1 = oldF1 - oldF2; + L F2 = - F1; + #else +* Otherwise, we must perform a full comparison. + + #ifdef `CacheWrite' +* Output the generated diagrams into a temporary file. + .sort:save-tmp; + Format 255; + #write <`cachetmpfile'> " %+E", F1 + Format 72; + #close <`cachetmpfile'> + #endif + + Skip; + + L F2 = + #switch `DiagramGenerator' + #case qgraf + #call RunQgraf(in=`in',out=`out',loops=`loops',options=`qgrafoptions') + #break + #case feyngraph + #call RunFeynGraph(in=`in',out=`out',loops=`loops',options=`feyngraphoptions') + #break + #endswitch + ; +#ifdef `Verbose' + P +sss; +#endif + .sort:run-`DiagramGenerator'; + #$nF2 = termsin_(F2); + + #call MakeDiagrams(`numIn',`numOut') + #call CompareDiagrams(F1,F2,F0) + + #ifdef `CacheWrite' +* If the generated diagrams are consistent with those by the other generator +* then write the cache file. + #if (`ZERO_F1' == 1) && (`ZERO_F2' == 1) + #message Write cache file: `cachefile' + #system mkdir -p "`CacheDir'" + #system echo "* Compared with `DiagramGeneratorVersion'" >"`cachetmpfile2'" + #system echo "* Number of diagrams: FORM: `$nF1', `DiagramGeneratorName': `$nF2'" >>"`cachetmpfile2'" + #system sed "s/^ */ /" "`cachetmpfile'" >>"`cachetmpfile2'" + #system mv "`cachetmpfile2'" "`cachefile'" + #endif + #remove <`cachetmpfile'> + #endif + #endif +#endprocedure + +*--#] DoComparison : +*--#] diagram_compare_include : +*--#[ qgraf_compare_include : +#- +#include compare.frm # diagram_compare_include + +*--#[ RunQgraf : + +** +* Runs QGRAF and reads its output. +* +* Example +* ------- +* #call SetupQcdModel() +* L F = +* #call RunQgraf(in=qua,out=qua,loops=1,options=nosnail) +* ; +* +* Example Output +* -------------- +* F = +* vx(QUA(-2),qua(1),glu(2))*vx(QUA(1),qua(-1),glu(2)); +* +#procedure RunQgraf(?a) + #define in "" + #define out "" + #define loops "" + #define loopmomentum "" + #define options "" + #define next "options" + + #do a = {`?a',} + #ifdef `a' + #if "`keepleft_(`a',3)'" == "in=" + #redefine next "in" + #redefine a "`takeleft_(`a',3)'" + #elseif "`keepleft_(`a',4)'" == "out=" + #redefine next "out" + #redefine a "`takeleft_(`a',4)'" + #elseif "`keepleft_(`a',6)'" == "loops=" + #redefine next "loops" + #redefine a "`takeleft_(`a',6)'" + #elseif "`keepleft_(`a',14)'" == "loop_momentum=" + #redefine next "loopmomentum" + #redefine a "`takeleft_(`a',14)'" + #elseif "`keepleft_(`a',8)'" == "options=" + #redefine next "options" + #redefine a "`takeleft_(`a',8)'" + #endif + #ifdef `a' + #ifdef ``next'' + #redefine `next' "``next'',`a'" + #else + #redefine `next' "`a'" + #endif + #endif + #endif + #enddo + + #ifndef `loopmomentum' + #redefine loopmomentum "k" + #endif + + #remove <`QgrafOutputFile'> + + #define filename "`QgrafInputFile'" + #create <`filename'> + #write <`filename'> "output = \'`QgrafOutputFile'\';" + #write <`filename'> "style = \'`QgrafStyleFile'\';" + #write <`filename'> "model = \'`QgrafModelFile'\';" + #write <`filename'> "in = `in';" + #write <`filename'> "out = `out';" + #write <`filename'> "loops = `loops';" + #write <`filename'> "loop_momentum = `loopmomentum';" + #write <`filename'> "options = `options';" + #close <`filename'> + + #system qgraf + +#ifdef `Verbose' + #system cat `QgrafOutputFile' +#endif + #include `QgrafOutputFile' +#endprocedure + +*--#] RunQgraf : +*--#[ SetupQgrafFormStyle : + +** +* Prepares the QGRAF style file for FORM. +* +#procedure SetupQgrafFormStyle() + #define filename "`QgrafStyleFile'" + #create <`filename'> + #write <`filename'> "" + #write <`filename'> "*" + #write <`filename'> "* file generated by " + #write <`filename'> "*" + #write <`filename'> "* *" + #write <`filename'> "*" + #write <`filename'> "#redefine QgrafVersion \"\"" + #write <`filename'> "*" + #write <`filename'> "" + #write <`filename'> "*--#[[ d:" + #write <`filename'> " " + #write <`filename'> " *vx((),)" + #write <`filename'> "" + #write <`filename'> "*--#]] d: " + #write <`filename'> "" + #write <`filename'> "*" + #write <`filename'> "* end" + #write <`filename'> "*" + #write <`filename'> "" + #close <`filename'> +#endprocedure + +*--#] SetupQgrafFormStyle : + +#call SetupQgrafFormStyle() +#redefine DiagramGenerator "qgraf" +#redefine DiagramGeneratorVersion "`~QgrafVersion'" + +*--#] qgraf_compare_include : +*--#[ feyngraph_compare_include : +#- +#include compare.frm # diagram_compare_include + +*--#[ RunFeynGraph : + +** +* Runs FeynGraph and reads its output. +* +* Example +* ------- +* #call SetupQcdModel() +* L F = +* #call RunFeynGraph(in=qua,out=qua,loops=1,options=self_loops=0) +* ; +* +* Example Output +* -------------- +* F = +* vx(QUA(-2),qua(3),glu(2))*vx(QUA(3),qua(-1),glu(2)); +* +#procedure RunFeynGraph(?a) + #define in "" + #define out "" + #define loops "" + #define loopmomentum "" + #define options "" + #define next "options" + + #do a = {`?a',} + #ifdef `a' + #if "`keepleft_(`a',3)'" == "in=" + #redefine next "in" + #redefine a "`takeleft_(`a',3)'" + #elseif "`keepleft_(`a',4)'" == "out=" + #redefine next "out" + #redefine a "`takeleft_(`a',4)'" + #elseif "`keepleft_(`a',6)'" == "loops=" + #redefine next "loops" + #redefine a "`takeleft_(`a',6)'" + #elseif "`keepleft_(`a',14)'" == "loop_momentum=" + #redefine next "loopmomentum" + #redefine a "`takeleft_(`a',14)'" + #elseif "`keepleft_(`a',8)'" == "options=" + #redefine next "options" + #redefine a "`takeleft_(`a',8)'" + #endif + #ifdef `a' + #ifdef ``next'' + #redefine `next' "``next'',`a'" + #else + #redefine `next' "`a'" + #endif + #endif + #endif + #enddo + + #ifndef `loopmomentum' + #redefine loopmomentum "k" + #endif + + #remove <`FeynGraphOutputFile'> + + #define filename "`FeynGraphInputFile'" + #create <`filename'> + #write <`filename'> "template = \"`FeynGraphTemplateFile'\"" + #write <`filename'> "[process]" + #write <`filename'> "in = [%" + #do a = {`in',} + #ifdef `a' + #write <`filename'> "\"`a'\", %" + #endif + #enddo + #write <`filename'> "]" + #write <`filename'> "out = [%" + #do a = {`out',} + #ifdef `a' + #write <`filename'> "\"`a'\", %" + #endif + #enddo + #write <`filename'> "]" + #write <`filename'> "loops = `loops'" + #write <`filename'> "model = \"`QgrafModelFile'\"" + #write <`filename'> "momenta = [%" + #define i "1" + #do a = {`in',} + #ifdef `a' + #write <`filename'> "\"p`i++'\", %" + #endif + #enddo + #redefine i "1" + #do a = {`out',} + #ifdef `a' + #write <`filename'> "\"q`i++'\", %" + #endif + #enddo + #do i = 1, `loops' + #write <`filename'> "\"k`i'\", %" + #enddo + #write <`filename'> "]" + #write <`filename'> "[filter]" + #do a = {`options',} + #ifdef `a' + #switch `a' + #case onshell + #write <`filename'> "`a' = true" + #break + #default + #write <`filename'> "`a'" + #endswitch + #endif + #enddo + #close <`filename'> + + #system feyngraph `FeynGraphInputFile' + #pipe echo \#redefine FeynGraphVersion \"$(feyngraph --version)\" + +#ifdef `Verbose' + #system cat `FeynGraphOutputFile' +#endif + #include `FeynGraphOutputFile' +#endprocedure + +*--#] RunFeynGraph : +*--#[ SetupFeynGraphTemplate : + +** +* Prepares the FeynGraph template file for FORM. +* +#procedure SetupFeynGraphTemplate() + #define filename "`FeynGraphTemplateFile'" + #create <`filename'> + #write <`filename'> "{\%- for diag in diagrams -\%}" + #write <`filename'> " {\%- set n_in = diag.n_in() -\%}" + #write <`filename'> " {\%- set n_ext = diag.n_ext() -\%}" + #write <`filename'> "*--#[[ d{{ loop.index }}:" + #write <`filename'> " {{ \"+\" ~ diag.sign() if diag.sign() >= 0 else diag.sign() -}}" + #write <`filename'> " {{ \"/\" ~ diag.symmetry_factor() if diag.symmetry_factor() != 1 else \"\" }}" + #write <`filename'> "{\%- for vertex in diag.vertices() \%}" + #write <`filename'> " *vx(" + #write <`filename'> "{\%- for par in vertex.particles_ordered() -\%}" + #write <`filename'> " {\%- set j = vertex.propagators_ordered()[loop.index0].id() -\%}" + #write <`filename'> " {\%- if j < n_in -\%}" + #write <`filename'> " {\%- set j = - 2 * j - 1 -\%}" + #write <`filename'> " {\%- elif j < n_ext -\%}" + #write <`filename'> " {\%- set j = - 2 * (j - n_in + 1) -\%}" + #write <`filename'> " {\%- endif -\%}" + #write <`filename'> " {{ par.name() }}({{ j }})" + #write <`filename'> " {\%- if not loop.last -\%},{\%- endif -\%}" + #write <`filename'> "{\%- endfor -\%}" + #write <`filename'> ")" + #write <`filename'> "{\%- endfor \%}" + #write <`filename'> "*--#]] d{{ loop.index }}: " + #write <`filename'> "{\% endfor \%}" + #close <`filename'> +#endprocedure + +*--#] SetupFeynGraphTemplate : + +#call SetupFeynGraphTemplate() +#redefine DiagramGenerator "feyngraph" +#redefine DiagramGeneratorVersion "feyngraph `~FeynGraphVersion'" + +*--#] feyngraph_compare_include : +*--#[ qgraf_phi4_phiphi_phiphi_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(phi4,in=phi,phi,out=phi,phi,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 7; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_phi4_phiphi_phiphi_1 : +*--#[ qgraf_phi4_phiphi_phiphi_2 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(phi4,in=phi,phi,out=phi,phi,loops=2,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 42; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_phi4_phiphi_phiphi_2 : +*--#[ qgraf_qcd_qua_qua_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=qua,out=qua,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 4; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_qua_qua_1 : +*--#[ qgraf_qcd_glu_glu_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=glu,out=glu,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 7; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_glu_glu_1 : +*--#[ qgraf_qcd_gho_gho_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=gho,out=gho,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 4; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_gho_gho_1 : +*--#[ qgraf_qcd_qua_quaglu_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=qua,out=qua,glu,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 17; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_qua_quaglu_1 : +*--#[ qgraf_qcd_glu_gluglu_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=glu,out=glu,glu,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 32; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_glu_gluglu_1 : +*--#[ qgraf_qcd_gluglu_gluglu_1 : +#include- compare.frm # qgraf_compare_include +#call DoComparison(qcd,in=glu,glu,out=glu,glu,loops=1,form_options=,qgraf_options=) +.end +assert succeeded? +assert nterms("F1", 0) == 223; +assert nterms("F1") == 0; +assert nterms("F2") == 0; +*--#] qgraf_qcd_gluglu_gluglu_1 :