From 75c613ccf6016c8fb3920bc43e71918c2eee5a49 Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Wed, 7 Jan 2026 10:52:44 +0000
Subject: [PATCH 1/9] fix(grcc): skip symmetry factors for topologiesonly_
 output

---
 sources/diawrap.cc | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/sources/diawrap.cc b/sources/diawrap.cc
index e9b0800f..ec981c4a 100644
--- a/sources/diawrap.cc
+++ b/sources/diawrap.cc
@@ -704,15 +704,6 @@ Bool ProcessTopology(EGraph *eg, void *ti)
 		*fill++ = 0; *fill++ = 1; *fill++ = 5;
 	}
 //
-//	Symmetry factors. We let Normalize do the multiplication.
-//
-	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;
-	}
-//
 //	finish it off
 //
 	while ( tail < tend ) *fill++ = *tail++;

From 59a743b7443dd682344d7fa3e119018dd235ffea Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Wed, 7 Jan 2026 11:33:53 +0000
Subject: [PATCH 2/9] fix(grcc): forbid numerical couplings in Vertex

The manual (currently) states that "The coupling constant is either
the number 1, or a (product of) symbol(s) to an integer power. If
the symbols have not been declared before they will be declared
during the reading of the vertex declaration.".

However, "1" causes empty arguments to be printed in node_ functions,
numbers other than "1" are accepted, and it doesn't work at all for
vertices with four or more legs.

Since the user can modify the couplings as they wish after diagram
generation with replace_ etc, demand symbols in the Vertex definition.
---
 sources/model.c | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/sources/model.c b/sources/model.c
index 81ce8606..45330963 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;

From 78ad49369a69ddf386f32d378fdfca64daaf3060 Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Wed, 7 Jan 2026 11:45:21 +0000
Subject: [PATCH 3/9] fix(grcc): forbid negative coupling powers

These are not accepted by grcc.cc anyway, so print a more useful
error and terminate.
---
 sources/model.c | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/sources/model.c b/sources/model.c
index 45330963..18026e0c 100644
--- a/sources/model.c
+++ b/sources/model.c
@@ -342,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;
 			}

From 46ace12b86da20f44a5f3fa5b33a3207d3217ead Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Thu, 8 Jan 2026 15:50:30 +0000
Subject: [PATCH 4/9] refactor(grcc): use GRCC_ABORT macro in all terminating
 cases

If compiling as part of FORM, call FORM's Terminate rather than exit.
---
 sources/grcc.cc     | 12 ++++++------
 sources/grccparam.h | 10 ++++++----
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/sources/grcc.cc b/sources/grcc.cc
index e2d2829a..a3e411a5 100644
--- a/sources/grcc.cc
+++ b/sources/grcc.cc
@@ -179,7 +179,7 @@ Options::Options(void)
         fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n");
         fprintf(GRCC_Stderr, "nOptDef=%d, GRCC_OPT_Size=%d\n",
                 nOptDef, GRCC_OPT_Size);
-        exit(1);
+        GRCC_ABORT();
     }
     
     model  = NULL;
@@ -201,7 +201,7 @@ Options::Options(void)
         fprintf(GRCC_Stderr, "*** Options: inconsistent default values\n");
         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++) {
@@ -6463,11 +6463,11 @@ void EGraph::fromDGraph(DGraph *dg)
 
     if (dg->nnodes > GRCC_MAXNODES) {
         fprintf(GRCC_Stderr, "*** too many nodes (GRCC_MAXNODES)\n");
-        exit(1);
+        GRCC_ABORT();
     }
     if (dg->nedges > GRCC_MAXEDGES) {
         fprintf(GRCC_Stderr, "*** too many edges (GRCC_MAXEDGES)\n");
-        exit(1);
+        GRCC_ABORT();
     }
 
     for (e = 0; e < dg->nedges; e++) {
@@ -6475,7 +6475,7 @@ void EGraph::fromDGraph(DGraph *dg)
             fprintf(GRCC_Stderr, "*** undefined node:");
             fprintf(GRCC_Stderr, "edge[%d] = {%d, %d}\n", 
                     e, dg->edges[e][0], dg->edges[e][0]);
-            exit(1);
+            GRCC_ABORT();
         }
     }
 
@@ -7041,7 +7041,7 @@ Bool EGraph::optQGrafM(Options *opt)
     }
     if (maxlegs >= GRCC_MAXEDGES) {
         printf("*** table overflow\n");
-        abort();
+        GRCC_ABORT();
     }
     for (int k = 0; k < GRCC_MAXEDGES; k++) {
         nopis[k] = 0;
diff --git a/sources/grccparam.h b/sources/grccparam.h
index 88d03e5d..b8333683 100644
--- a/sources/grccparam.h
+++ b/sources/grccparam.h
@@ -31,10 +31,12 @@
 #define GRCC_Stderr          stderr
 */
 #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

From f48061b5685a05654d26d3d814eabe4e803249ee Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Thu, 8 Jan 2026 16:37:06 +0000
Subject: [PATCH 5/9] refactor(grcc): wrap prints to stdout and stderr

When compiling under FORM, we can format the output with vsnprintf,
and then just give the formatted string to MesPrint, such that it ends
up in the FORM log file.

When compiling without FORM, just forward the format string and the
args to vfprintf, which prints them as previously.
---
 sources/grcc.cc     | 912 +++++++++++++++++++++++---------------------
 sources/grccparam.h |   6 +-
 2 files changed, 474 insertions(+), 444 deletions(-)

diff --git a/sources/grcc.cc b/sources/grcc.cc
index a3e411a5..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,8 +177,8 @@ 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);
         GRCC_ABORT();
     }
@@ -198,8 +199,8 @@ 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);
         GRCC_ABORT();
     }
@@ -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,18 +6463,18 @@ void EGraph::fromDGraph(DGraph *dg)
     int n0, n1, e;
 
     if (dg->nnodes > GRCC_MAXNODES) {
-        fprintf(GRCC_Stderr, "*** too many nodes (GRCC_MAXNODES)\n");
+        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");
+        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]);
             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,7 +7041,7 @@ Bool EGraph::optQGrafM(Options *opt)
         maxlegs = Max(maxlegs, econn->opics[j].nlegs);
     }
     if (maxlegs >= GRCC_MAXEDGES) {
-        printf("*** table overflow\n");
+        grcc_fprintf(GRCC_Stdout, "*** table overflow\n");
         GRCC_ABORT();
     }
     for (int k = 0; k < GRCC_MAXEDGES; k++) {
@@ -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 b8333683..7ab17eb2 100644
--- a/sources/grccparam.h
+++ b/sources/grccparam.h
@@ -26,10 +26,8 @@
 
 #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
 
 #ifndef NOFORM

From d6eb8fa4f490457c5bd57188b39ea0040561e8c5 Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Mon, 12 Jan 2026 16:32:47 +0000
Subject: [PATCH 6/9] fix: allow empty momenta sets in diagrams_

This means we don't have to provide "dummy" external momenta when
generating vacuum graphs. An empty internal momentum set is more
dubious, but allow it in the input nonetheless.
---
 sources/proces.c | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/sources/proces.c b/sources/proces.c
index 40fb5cf5..401f9f7a 100644
--- a/sources/proces.c
+++ b/sources/proces.c
@@ -1220,8 +1220,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

From 3b4cc630320c896f15777dc587c2eb115398f9fa Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Tue, 13 Jan 2026 11:46:30 +0000
Subject: [PATCH 7/9] fix: verify sufficient internal and external momenta
 passed to diagrams_

When generating diagrams, ensure that sufficient momenta have been
provided for external and internal lines, to avoid later crashes
during set element lookup.
---
 sources/diawrap.cc | 37 +++++++++++++++++++++++++++++++++++--
 1 file changed, 35 insertions(+), 2 deletions(-)

diff --git a/sources/diawrap.cc b/sources/diawrap.cc
index ec981c4a..360b1256 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,6 +721,15 @@ Bool ProcessTopology(EGraph *eg, void *ti)
 		*fill++ = 0; *fill++ = 1; *fill++ = 5;
 	}
 //
+//	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++;
@@ -838,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++ ) {

From 3168a92c37906d99190a7398dd4ec7f981d1ada8 Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Mon, 12 Jan 2026 10:40:04 +0000
Subject: [PATCH 8/9] remove the topologies_ function and associated code

Clean up code related to topologies_, now that we only use the
TopologiesOnly_ option of diagrams_. Clean up old grcc code which
generates topologies only.
---
 sources/declare.h   |    2 -
 sources/diagrams.c  |   71 --
 sources/diawrap.cc  |  110 ---
 sources/gentopo.cc  | 1583 -------------------------------------------
 sources/gentopo.h   |  162 -----
 sources/proces.c    |   32 +-
 sources/topowrap.cc |  283 --------
 7 files changed, 5 insertions(+), 2238 deletions(-)
 delete mode 100644 sources/gentopo.cc
 delete mode 100644 sources/gentopo.h
 delete mode 100644 sources/topowrap.cc

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 360b1256..458c6617 100644
--- a/sources/diawrap.cc
+++ b/sources/diawrap.cc
@@ -1000,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 <config.h>
-#endif
-#endif
-
-#include <iostream>
-#include <iomanip>
-#include <sstream>
-#include <cstdio>
-#include <stdlib.h>
-
-#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/proces.c b/sources/proces.c
index 401f9f7a..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 ) {
 /*
@@ -1262,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;
@@ -1288,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;
@@ -4080,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 <http://www.gnu.org/licenses/>.
- */
-/*
-  	#] 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 : 
-*/
-

From dde9a81a80612cde916cd6c5fef5761251fb65da Mon Sep 17 00:00:00 2001
From: Josh Davies <jodavies@users.noreply.github.com>
Date: Tue, 13 Jan 2026 12:41:03 +0000
Subject: [PATCH 9/9] test: add cases for new diagrams_ related error messages

---
 check/features.frm | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

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 :