===================================================================
RCS file: /home/cvs/OpenXM/src/phc/phc6.c,v
retrieving revision 1.1.1.1
retrieving revision 1.7
diff -u -p -r1.1.1.1 -r1.7
--- OpenXM/src/phc/phc6.c	1999/10/08 02:12:14	1.1.1.1
+++ OpenXM/src/phc/phc6.c	2019/01/01 11:06:10	1.7
@@ -1,4 +1,5 @@
 /*  phc6.c ,  yama:1999/sm1-prog/phc6.c */
+/* $OpenXM: OpenXM/src/phc/phc6.c,v 1.6 2019/01/01 07:17:52 takayama Exp $ */
 /* This is a simple C interface to the black-box solver of phc.
 ** Requirements:
 **  1) executable version of phc will be searched in the following order:
@@ -13,6 +14,8 @@
 #include <sys/stat.h>
 #include <unistd.h>
 #include <stdlib.h>
+#include <string.h>
+#include  "gc.h"
 
 /* Definition of class identifiers. */
 #define Snull             0
@@ -25,7 +28,7 @@ union cell {
   int ival;
   char *str;
   struct phc_object *op;
-  long double longdouble;
+  double longdouble;
 };
 struct phc_object{
   int tag;                /* class identifier */
@@ -36,17 +39,18 @@ struct phc_object{
 /* Memory allocation function.  
    Use your favorite memory allocation function.
    I recommend not to use malloc and to use gc4.14 for large applications. */
-#define sGC_malloc(n)   malloc(n)
+//#define sGC_malloc(n)   malloc(n)
+#define sGC_malloc(n)   GC_malloc(n)
 
 /********** macros to use Sarray **************/
 /* put to Object Array */
 #define phc_putoa(ob,i,cc) {\
-if ((ob).tag != Sarray) {fprintf(stderr,"Warning: PUTOA is for an array of objects\n");} else \
-{if ((0 <= (i)) && ((i) < (ob).lc.ival)) {\
-  (ob.rc.op)[i] = cc;\
-}else{\
-  fprintf(stderr,"Warning: PUTOA, the size is %d.\n",(ob).lc.ival);\
-}}}
+                   if ((ob).tag != Sarray) {fprintf(stderr,"Warning: PUTOA is for an array of objects\n");} else \
+                     {if ((0 <= (i)) && ((i) < (ob).lc.ival)) {\
+                      (ob.rc.op)[i] = cc;\
+                    }else{\
+                       fprintf(stderr,"Warning: PUTOA, the size is %d.\n",(ob).lc.ival);\
+                  }}}
 #define phc_getoa(ob,i) ((ob.rc.op)[i])
 #define phc_getoaSize(ob) ((ob).lc.ival)
 
@@ -55,38 +59,67 @@ struct phc_object phc_newObjectArray(int size);
 void phc_printObject(FILE *fp,struct phc_object ob);
 char *phc_generateUniqueFileName(char *s);
 char *phc_which(char *s);   /* search a path for the file s */
-struct phc_object phc_complexTo(long double r, long double i);
+struct phc_object phc_complexTo(double r, double i);
 
 
-int phc_scan_for_string(FILE *fp, char str[], int lenstr);
+int phc_scan_for_string(FILE *fp, char str[]);
 struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim );
-struct phc_object phc_scan_output_of_phc(char *fname);
+struct phc_object phc_scan_output_of_phc(FILE *fp);
 struct phc_object phc_call_phc(char *sys); 
 
 int phc_verbose = 0;
 int phc_overwrite = 1;  /* Always use tmp.input.0 and tmp.output.0
-			   for work files. */
+                           for work files. */
 
-main(int argc, char *argv[]) {
+int main(int argc, char *argv[]) {
   struct phc_object ob;
+  struct phc_object ob2;
+  int n2;
   int n,i,dim;
 #define INPUTSIZE 4096
   char input[INPUTSIZE];
+  char fname[INPUTSIZE];
 #define A_SIZE 1024
   char a[A_SIZE];
   int message = 0;
+  FILE *fp;
   for (i=1; i<argc; i++) {
     if (strcmp(argv[i],"-v") == 0) {
-      phc_verbose = 1;
+      phc_verbose = 1; message=1;
     }else if (strcmp(argv[i],"-g") == 0) {
       phc_overwrite = 0;
+    }else if (strcmp(argv[i],"-file") == 0) {
+      /* For debugging. */
+      i++;
+      strncpy(fname,argv[i],INPUTSIZE-1);
+      fp = fopen(fname,"r");
+      if (fp == NULL) {
+	fprintf(stderr,"File %s is not found.\n",fname);
+        return(-1);
+      }
+      ob = phc_scan_output_of_phc(fp);  /* stable mixed vol,  out of torus */
+      n = phc_getoaSize(ob);
+      ob2 = phc_scan_output_of_phc(fp); /* in torus */
+      n2 = phc_getoaSize(ob2);
+      if (phc_verbose) fprintf(stderr,"n=%d, n2=%d\n",n,n2);
+      printf("[\n");
+      for (i=0; i<n; i++) {
+        phc_printObject(stdout,phc_getoa(ob,i));
+        if ((i != n-1) || (n2 > 0)) printf(" ,\n"); else printf(" \n");
+      }
+      for (i=0; i<n2; i++) {
+        phc_printObject(stdout,phc_getoa(ob2,i));
+        if (i != n2-1) printf(" ,\n"); else printf(" \n");
+      }
+      printf("]\n");
+      return(0);
     }else if (strcmp(argv[i],"-i") == 0) {
       ob = phc_call_phc(argv[i+1]);
       n = phc_getoaSize(ob);
       printf("[\n");
       for (i=0; i<n; i++) {
-	phc_printObject(stdout,phc_getoa(ob,i));
-	if (i != n-1) printf(" ,\n"); else printf(" \n");
+        phc_printObject(stdout,phc_getoa(ob,i));
+        if (i != n-1) printf(" ,\n"); else printf(" \n");
       }
       printf("]\n");
       exit(0);
@@ -98,21 +131,22 @@ main(int argc, char *argv[]) {
   }
   while (1) {
     if (message) printf("dim= ");
-    if (scanf("%d",&dim)<0) break;
+    if (fgets(input,INPUTSIZE,stdin) <= 0) break;
+    sscanf(input,"%d",&dim);
     sprintf(input,"%d\n",dim);
     if (message) printf("Input %d equations please.\n",dim);
     for (i=0; i<dim; i++) {
-     if (message) {printf("eq[%d] = ",i);  fflush(stdout);}
-     do {
-       fgets(a,A_SIZE-1, stdin);
-     } while (strlen(a) == 0);
-     if (strlen(a) >= A_SIZE-3) {
-       fprintf(stderr,"Too big input for the input buffer a.\n"); exit(10);
-     }
-     if (strlen(input)+strlen(a) >= INPUTSIZE) {
-       fprintf(stderr,"Too big input for the input buffer input.\n"); exit(10);
-     }
-     sprintf(input+strlen(input),"%s\n",a);
+      if (message) {printf("eq[%d] = ",i);  fflush(stdout);}
+      do {
+        fgets(a,A_SIZE-1, stdin);
+      } while (strlen(a) == 0);
+      if (strlen(a) >= A_SIZE-3) {
+        fprintf(stderr,"Too big input for the input buffer a.\n"); exit(10);
+      }
+      if (strlen(input)+strlen(a) >= INPUTSIZE) {
+        fprintf(stderr,"Too big input for the input buffer input.\n"); exit(10);
+      }
+      sprintf(input+strlen(input),"%s\n",a);
     }
     ob = phc_call_phc(input);
     if (message) {
@@ -126,52 +160,91 @@ main(int argc, char *argv[]) {
     }
     printf("]\n");
   }
+  return(0);
 }
 
-int phc_scan_for_string(FILE *fp, char str[], int lenstr)
+int phc_scan_for_string(FILE *fp, char s[])
      /*
   **  Scans the file fp for a certain string str of length lenstr+1.
   **  Reading stops when the string has been found, then the variable
   **  on return equals 1, otherwise 0 is returned.
   */
+#define BUF_SIZE 1024
 {
-  char buf[lenstr+1];
+  char buf[BUF_SIZE];
+  int pt,n,c,i;
+  pt = 0;
+  n=strlen(s);
+  if (n > BUF_SIZE-2) {
+	fprintf(stderr,"Too long string for scan_for_string.\n");
+	exit(1);
+  }
+  for (i=0; i<n; i++) {
+	buf[i] = fgetc(fp); buf[i+1]='\0';
+	if (buf[i] == EOF) return 0;
+  }
+  if (strcmp(s,buf) == 0) return 1;
+  while ((c=fgetc(fp)) != EOF) {
+	for (i=0; i<n; i++) {
+	  buf[i] = buf[i+1];
+	}
+	buf[n-1]=c; buf[n] = '\0';
+	if (strcmp(s,buf) == 0) return 1;
+  }
+  return 0;
+}
+
+int phc_scan_for_string_old(FILE *fp, char str[])
+     /*
+  **  Scans the file fp for a certain string str of length lenstr+1.
+  **  Reading stops when the string has been found, then the variable
+  **  on return equals 1, otherwise 0 is returned.
+  */
+{
+#define BUF_SIZE 1024
+  char buf[BUF_SIZE];
   char ch;
   int index,i,compare,npaths,dim,found;
+  int lenstr;
+  lenstr = strlen(str);
+  if (lenstr >= BUF_SIZE-1) {
+    fprintf(stderr,"Too long string in phc_scan_for_string\n");
+    exit(-1);
+  }
   index = -1;
   found = 0;
-  while ((fscanf(fp,"%c",&ch)!=EOF) && found == 0)
+  while (((ch = fgetc(fp))!=EOF) && found == 0)
     {
       if (index == -1 && ch == str[0])
-	{
-	  index = 0;
-	  buf[index] = ch;
-	}
+        {
+          index = 0;
+          buf[index] = ch;
+        }
       else
-	{
-	  if (index == lenstr)
-	    {
-	      compare = 0;
-	      for (i=0; i<lenstr+1; i++)
-		{
-		  if (buf[i]!=str[i])
-		    {
-		      compare = compare+1;
-		    }
-		}
-	      if (compare == 0)
-		{
-		  found = 1;
-		}
-	      index = -1;
-	    }
-	  else
-	    if (index > -1 && index < lenstr)
-	      {
-		index = index+1;
-		buf[index] = ch;
-	      }
-	}
+        {
+          if (index == lenstr)
+            {
+              compare = 0;
+              for (i=0; str[i] != '\0'; i++)
+                {
+                  if (buf[i]!=str[i])
+                    {
+                      compare = compare+1;
+                    }
+                }
+              if (compare == 0)
+                {
+                  found = 1;
+                }
+              index = -1;
+            }
+          else
+            if (index > -1 && index < lenstr)
+              {
+                index = index+1;
+                buf[index] = ch;
+              }
+        }
       if (found == 1) break;
     }
   return found;
@@ -183,57 +256,92 @@ struct phc_object phc_scan_solutions(FILE *fp, int npa
   **  The tolerance for the residual to a solution is set to 1.0E-12.
   **  Returns solutions.
   */
+#define BUFSIZE 1024
 {
   struct phc_object rob,sob;
   char ch;
   int fnd,i,j,nsols;
-  float res;
-  long double realpart;
-  long double imagpart;
-  long double realparts[npaths][dim];
-  long double imagparts[npaths][dim];
+  double res;
+  double realpart;
+  double imagpart;
+  /*  double realparts[npaths][dim];
+  double imagparts[npaths][dim]; */
+  double *realparts;
+  double *imagparts;
+  double *res_list;
+  char buf[BUFSIZE];
   nsols = 0;
-  while (fscanf(fp,"%c",&ch)!=EOF) 
+  realparts = (double *)sGC_malloc(sizeof(double)*(npaths*dim+1));
+  imagparts = (double *)sGC_malloc(sizeof(double)*(npaths*dim+1));
+  res_list = (double *)sGC_malloc(sizeof(double)*npaths+1);
+  while ((ch = fgetc(fp)) != EOF) 
     {
-      fnd = phc_scan_for_string(fp,"start residual :",15);
+      fnd = phc_scan_for_string(fp,"start residual :");
       if (fnd==1)
-	{
-	  fscanf(fp,"%E",&res);
-	  /* printf(" residual = "); printf("%E\n",res); */
-	  if (res < 1.0E-12) nsols = nsols+1;
-	  fnd = phc_scan_for_string(fp,"the solution for t :",19);
-	  for (i=0;i<dim;i++)
-	    {
-	      fnd = phc_scan_for_string(fp,":",0);
-	      fscanf(fp,"%LE",&realpart);
-	      fscanf(fp,"%LE",&imagpart);
-	      if (res < 1.0E-12)
-		{
-		  realparts[nsols-1][i] = realpart;
-		  imagparts[nsols-1][i] = imagpart;
-		}
-	    }
-	}
+        {
+	  fgets(buf,BUFSIZE-1,fp);
+          sscanf(buf,"%lg",&res);
+	  if (phc_verbose) {
+	    fprintf(stderr,"res in string: %s\n",buf);
+	    fprintf(stderr," residual = "); fprintf(stderr,"%le\n",res);
+	  }
+          nsols++;
+          /*
+          if (res < 1.0E-12) {
+            nsols = nsols+1;
+            if (nsols > npaths) {
+              fprintf(stderr,"Something is wrong in phc_scan_solutions\n");
+              fprintf(stderr,"npaths=%d, nsols=%d \n",npaths,nsols);
+              exit(-1);
+            }
+          }
+          */
+          fnd = phc_scan_for_string(fp,"the solution for t :");
+          res_list[nsols-1] = res;
+          for (i=0;i<dim;i++)
+            {
+              fnd = phc_scan_for_string(fp,":");
+	      fgets(buf,BUFSIZE-1,fp);
+	      sscanf(buf,"%lf %lf",&realpart,&imagpart);
+	      if (phc_verbose) fprintf(stderr,"sols in string: %s\n",buf);
+              /*fscanf(fp,"%le",&realpart);
+		fscanf(fp,"%le",&imagpart);*/
+              realparts[(nsols-1)*dim+i] = realpart;
+              imagparts[(nsols-1)*dim+i] = imagpart;
+            }
+	  if (nsols >= npaths) break;
+        }
     }
+  if (nsols < npaths) {
+    fprintf(stderr,"Warning: nsols < npaths. Check tmp.output.*\n");
+  }
   if(phc_verbose) fprintf(stderr,"  number of solutions = %i\n",nsols); 
   rob = phc_newObjectArray(nsols);
   for (i=0;i<nsols;i++)
     {
-      /* fprintf(stderr,"Solution %i :\n",i+1); */
+      if (phc_verbose) fprintf(stderr,"Solution %i :\n",i+1);
       sob = phc_newObjectArray(dim);
       for (j=0;j<dim;j++)
-	{
-	  /*
-	  printf("%20.14LE",realparts[i][j]); printf("  ");
-	  printf("%20.14LE",imagparts[i][j]); printf("\n");
-	  */
-	  phc_putoa(sob,j,phc_complexTo(realparts[i][j],imagparts[i][j]));
-	}
+        {
+          /*
+            printf("%20.14le",realparts[i][j]); printf("  ");
+            printf("%20.14le",imagparts[i][j]); printf("\n");
+            */
+		  if (phc_verbose) {
+            printf("%20.14le",realparts[i*dim+j]); printf("  ");
+            printf("%20.14le",imagparts[i*dim+j]); printf("\n");
+          }
+
+          /*phc_putoa(sob,j,phc_complexTo(realparts[i][j],imagparts[i][j]));*/
+          phc_putoa(sob,j,phc_complexTo(realparts[i*dim+j],
+                                        imagparts[i*dim+j]));
+        }
       phc_putoa(rob,i,sob);
     }
+  /* Todo, res_list[] is not returned. */
   return(rob);
 }
-struct phc_object phc_scan_output_of_phc(char *fname)
+struct phc_object phc_scan_output_of_phc(FILE *fp)
      /*
   **  Scans the file "output" of phc in two stages to get
   **   1) the number of paths and the dimension;
@@ -243,9 +351,18 @@ struct phc_object phc_scan_output_of_phc(char *fname)
   FILE *otp;
   char ch;
   int fnd,npaths,dim,i,nsols;
-  otp = fopen(fname,"r");
-  if (phc_verbose) fprintf(stderr,"Scanning the %s of phc.\n",fname);
-  fnd = phc_scan_for_string(otp,"THE SOLUTIONS :",14);
+  struct phc_object rob;
+  otp = fp;
+  if (otp == NULL) {
+    fprintf(stderr,"Error: the file is not found.\n");
+    exit(0);
+  }
+  fnd = phc_scan_for_string(otp,"THE SOLUTIONS :");
+  if (fnd < 1) {
+    if (phc_verbose) fprintf(stderr,"No more THE SOLUTIONS\n");
+    rob = phc_newObjectArray(0);
+    return(rob);
+  }
   fscanf(otp,"%i",&npaths);
   if (phc_verbose) fprintf(stderr,"  number of paths traced = %i\n",npaths);
   fscanf(otp,"%i",&dim);
@@ -260,6 +377,11 @@ struct phc_object phc_call_phc(char *sys)  /* call phc
   char *w;
   struct phc_object phc_NullObject ;
   struct stat statbuf;
+  FILE *fp;
+  struct phc_object ob;
+  struct phc_object ob2;
+  struct phc_object ob_all;
+  int i,n,n2;
 
   phc_NullObject.tag = Snull;
   f = phc_generateUniqueFileName("tmp.input");
@@ -270,7 +392,7 @@ struct phc_object phc_call_phc(char *sys)  /* call phc
     system(cmd);
   }
   inp = fopen(f,"w");
-  fprintf(inp,sys);
+  fprintf(inp,"%s",sys);
   fclose(inp);
   if ((w = phc_which("phc")) != NULL) {
     sprintf(cmd,"%s -b %s %s",w,f,outf);
@@ -284,13 +406,22 @@ struct phc_object phc_call_phc(char *sys)  /* call phc
     return(phc_NullObject);
   }else{
     if (phc_verbose) fprintf(stderr,"See the file %s for results.\n",outf);
-    return(phc_scan_output_of_phc(outf));
+    fp = fopen(outf,"r");
+    ob=phc_scan_output_of_phc(fp);
+    n = phc_getoaSize(ob);
+    ob2=phc_scan_output_of_phc(fp);
+    n2 = phc_getoaSize(ob2);
+    ob_all=phc_newObjectArray(n+n2);
+    if (phc_verbose) fprintf(stderr,"n=%d, n2=%d\n",n,n2);
+    for (i=0; i<n; i++) phc_putoa(ob_all,i,phc_getoa(ob,i));
+    for (i=0; i<n2; i++) phc_putoa(ob_all,n+i,phc_getoa(ob2,i));
+    return(ob_all);
   }
 }
 
 
 struct phc_object phc_newObjectArray(size) 
-int size;
+     int size;
 {
   struct phc_object rob;
   struct phc_object *op;
@@ -323,8 +454,8 @@ void phc_printObject(FILE *fp,struct phc_object ob)
     fprintf(fp,"]");
   }else if (ob.tag == SlongdoubleComplex) {
     /* Try your favorite way to print complex numbers. */
-    /*fprintf(fp,"(%20.14LE)+I*(%20.14LE)",ob.lc.longdouble,ob.rc.longdouble);*/
-    fprintf(fp,"[%Lf , %Lf]",ob.lc.longdouble,ob.rc.longdouble);
+    /*fprintf(fp,"(%20.14le)+I*(%20.14le)",ob.lc.longdouble,ob.rc.longdouble);*/
+    fprintf(fp,"[%.16lg , %.16lg]",ob.lc.longdouble,ob.rc.longdouble);
   }else{
     fprintf(stderr,"Unknown phc_object tag %d",ob.tag);
   }
@@ -355,7 +486,7 @@ char *phc_which(char *s) {
   a = getenv("OpenXM_HOME");
   if (a != NULL) {
     cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)+strlen(a)
-					    +strlen("/usr/local/bin/")+1));
+                                            +strlen("/usr/local/bin/")+1));
     if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
     strcpy(cmd,a); strcat(cmd,"/bin/"); strcat(cmd,s);
     if (stat(cmd,&statbuf) == 0) {
@@ -363,7 +494,7 @@ char *phc_which(char *s) {
     }
   }
   cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)
-					  +strlen("/usr/local/bin/")+1));
+                                          +strlen("/usr/local/bin/")+1));
   if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
   strcpy(cmd,"/usr/local/bin/"); strcat(cmd,s);
   if (stat(cmd,&statbuf) == 0) {
@@ -377,7 +508,7 @@ char *phc_which(char *s) {
 }
 
    
-struct phc_object phc_complexTo(long double r, long double i)
+struct phc_object phc_complexTo(double r, double i)
 {
   struct phc_object rob;
   rob.tag = SlongdoubleComplex;