Grcc fixes

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 :

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);

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

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 :