/* This is a simple C interface to the black-box solver of phc,
** written by Nobuki Takayama.
** Requirements:
**  1) executable version of phc must be on /tmp;
**  2) user of this program has write permissions to create
**     the files "input" and "output" in the directory where
**     this program is executed.
*/

#include <stdio.h>
#include <sys/stat.h>
#include <unistd.h>

/* Definition of class identifiers. */
#define SstringObject     5
#define Sarray            6  

/* Definition of Object */
union cell {
  int ival;
  char *str;
  struct object *op;
};
struct object{
  int tag;                /* class identifier */
  union cell lc;          /* left cell */
  union cell rc;          /* right cell */
};

/* Memory allocation function.  
   I recommend not to use malloc and to use gc4.14 for large applications. */
#define sGC_malloc(n)   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);\
}}}
#define phc_getoa(ob,i) ((ob.rc.op)[i])
#define phc_getoaSize(ob) ((ob).lc.ival)


/* prototypes */
struct object phc_newObjectArray(int size);
void phc_printObject(FILE *fp,struct object ob);
char *phc_generateUniqueFileName(char *s);
struct object phc_complexTo(long double r, long double i);
struct object phc_longdoubleToStringObject(long double r);

int phc_scan_for_string(FILE *fp, char str[], int lenstr);
struct object phc_scan_solutions(FILE *fp, int npaths, int dim );
struct object phc_scan_output_of_phc(char *fname);
struct object phc_call_phc(char *sys); 

main() {
  struct object ob;
  ob = phc_call_phc("2\n x**2 + y**2 - 1;\n x**2 + y**2 - 8*x - 3;\n");
  printf("-----------------------------------------------------------\n");
  phc_printObject(stdout,ob);
  printf("\n");
}

int phc_scan_for_string(FILE *fp, char str[], int lenstr)
     /*
  **  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.
  */
{
  char buf[lenstr+1];
  char ch;
  int index,i,compare,npaths,dim,found;
  index = -1;
  found = 0;
  while ((fscanf(fp,"%c",&ch)!=EOF) && found == 0)
    {
      if (index == -1 && ch == str[0])
	{
	  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 (found == 1) break;
    }
  return found;
}
struct object phc_scan_solutions(FILE *fp, int npaths, int dim )
     /*
  **  Scans the file for the solutions, from a list of length npaths,
  **  of complex vectors with dim entries.
  **  The tolerance for the residual to a solution is set to 1.0E-12.
  **  Returns the number of solutions.
  */
{
  struct 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];
  nsols = 0;
  while (fscanf(fp,"%c",&ch)!=EOF) 
    {
      fnd = phc_scan_for_string(fp,"start residual :",15);
      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;
		}
	    }
	}
    }
  /* 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); */
      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]));
	}
      phc_putoa(rob,i,sob);
    }
  return(rob);
}
struct object phc_scan_output_of_phc(char *fname)
     /*
  **  Scans the file "output" of phc in two stages to get
  **   1) the number of paths and the dimension;
  **   2) the solutions, vectors with residuals < 1.0E-12.
  */
{
  FILE *otp;
  char ch;
  int fnd,npaths,dim,i,nsols;
  otp = fopen(fname,"r");
  fprintf(stderr,"Scanning the %s of phc.\n",fname);
  fnd = phc_scan_for_string(otp,"THE SOLUTIONS :",14);
  fscanf(otp,"%i",&npaths);
  fprintf(stderr,"  number of paths traced = %i\n",npaths);
  fscanf(otp,"%i",&dim);
  fprintf(stderr,"  dimension of solutions = %i\n",dim);
  return(phc_scan_solutions(otp,npaths,dim));
}
struct object phc_call_phc(char *sys)  /* call phc, system as string */
{
  FILE *inp;
  char *f,*outf;
  char cmd[1024];
  f = phc_generateUniqueFileName("tmp.input");
  outf = phc_generateUniqueFileName("tmp.output");
  fprintf(stderr,"Creating file with name %s.\n",f);
  inp = fopen(f,"w");
  fprintf(inp,sys);
  fclose(inp);
  sprintf(cmd,"/tmp/phc -b %s %s",f,outf);
  fprintf(stderr,"Calling %s, black-box solver of phc.\n",cmd);
  system(cmd);
  fprintf(stderr,"See the file %s for results.\n",outf);
  return(phc_scan_output_of_phc(outf));
}


struct object phc_newObjectArray(size) 
int size;
{
  struct object rob;
  struct object *op;
  if (size > 0) {
    op = (struct object *)sGC_malloc(size*sizeof(struct object));
    if (op == (struct object *)NULL) {fprintf(stderr,"No memory\n");exit(1);}
  }else{
    op = (struct object *)NULL;
  }
  rob.tag = Sarray;
  rob.lc.ival = size;
  rob.rc.op = op;
  return(rob);
}

void phc_printObject(FILE *fp,struct object ob)
{
  int n,i;
  if (ob.tag == SstringObject) {
    fprintf(fp,"%s",ob.lc.str);
  }else if (ob.tag == Sarray) {
    n = phc_getoaSize(ob);
    fprintf(fp,"[");
    for (i=0; i<n; i++) {
      phc_printObject(fp,phc_getoa(ob,i));
      if (i <n-1) fprintf(fp,", ");
    }
    fprintf(fp,"]");
  }else{
    fprintf(stderr,"Unknown object tag %d",ob.tag);
  }
}
    

char *phc_generateUniqueFileName(char *s)
{
  char *t;
  int i;
  struct stat statbuf;
  t = (char *)sGC_malloc(sizeof(char)*strlen(s)+4+2);
  for (i=0; i<1000; i++) {
    /* Give up if we failed for 1000 names. */
    sprintf(t,"%s.%d",s,i);
    if (stat(t,&statbuf) < 0) {
      return(t);
    }
  }
  fprintf(stderr,"Could not generate a unique file name.");
  return(NULL);
}

struct object phc_longdoubleToStringObject(long double r)
{
  struct object rob;
  rob.tag = SstringObject;
  rob.lc.str = (char *) sGC_malloc(sizeof(char)*50); /*??*/
  sprintf(rob.lc.str,"%20.14LE",r);
  return(rob);
}
   
struct object phc_complexTo(long double r, long double i)
{
  struct object rob;
  rob = phc_newObjectArray(2);
  phc_putoa(rob,0,phc_longdoubleToStringObject(r));
  phc_putoa(rob,1,phc_longdoubleToStringObject(i));
  return(rob);
}