Annotation of OpenXM/src/phc/phc6.c, Revision 1.1
1.1 ! maekawa 1: /* phc6.c , yama:1999/sm1-prog/phc6.c */
! 2: /* This is a simple C interface to the black-box solver of phc.
! 3: ** Requirements:
! 4: ** 1) executable version of phc will be searched in the following order:
! 5: ** OpenXM_HOME/bin/, /usr/local/bin/phc, /tmp/phc, your command search path.
! 6: ** Here, PHC_LIBDIR is an environment variable.
! 7: ** 2) user of this program has write permissions to create
! 8: ** the files "tmp.input.xx" and "tmp.output.xx" in the directory where
! 9: ** this program is executed. xx are numbers.
! 10: */
! 11:
! 12: #include <stdio.h>
! 13: #include <sys/stat.h>
! 14: #include <unistd.h>
! 15: #include <stdlib.h>
! 16:
! 17: /* Definition of class identifiers. */
! 18: #define Snull 0
! 19: #define SstringObject 5
! 20: #define Sarray 6
! 21: #define SlongdoubleComplex 601
! 22:
! 23: /* Definition of Object */
! 24: union cell {
! 25: int ival;
! 26: char *str;
! 27: struct phc_object *op;
! 28: long double longdouble;
! 29: };
! 30: struct phc_object{
! 31: int tag; /* class identifier */
! 32: union cell lc; /* left cell */
! 33: union cell rc; /* right cell */
! 34: };
! 35:
! 36: /* Memory allocation function.
! 37: Use your favorite memory allocation function.
! 38: I recommend not to use malloc and to use gc4.14 for large applications. */
! 39: #define sGC_malloc(n) malloc(n)
! 40:
! 41: /********** macros to use Sarray **************/
! 42: /* put to Object Array */
! 43: #define phc_putoa(ob,i,cc) {\
! 44: if ((ob).tag != Sarray) {fprintf(stderr,"Warning: PUTOA is for an array of objects\n");} else \
! 45: {if ((0 <= (i)) && ((i) < (ob).lc.ival)) {\
! 46: (ob.rc.op)[i] = cc;\
! 47: }else{\
! 48: fprintf(stderr,"Warning: PUTOA, the size is %d.\n",(ob).lc.ival);\
! 49: }}}
! 50: #define phc_getoa(ob,i) ((ob.rc.op)[i])
! 51: #define phc_getoaSize(ob) ((ob).lc.ival)
! 52:
! 53: /* prototypes */
! 54: struct phc_object phc_newObjectArray(int size);
! 55: void phc_printObject(FILE *fp,struct phc_object ob);
! 56: char *phc_generateUniqueFileName(char *s);
! 57: char *phc_which(char *s); /* search a path for the file s */
! 58: struct phc_object phc_complexTo(long double r, long double i);
! 59:
! 60:
! 61: int phc_scan_for_string(FILE *fp, char str[], int lenstr);
! 62: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim );
! 63: struct phc_object phc_scan_output_of_phc(char *fname);
! 64: struct phc_object phc_call_phc(char *sys);
! 65:
! 66: int phc_verbose = 0;
! 67: int phc_overwrite = 1; /* Always use tmp.input.0 and tmp.output.0
! 68: for work files. */
! 69:
! 70: main(int argc, char *argv[]) {
! 71: struct phc_object ob;
! 72: int n,i,dim;
! 73: #define INPUTSIZE 4096
! 74: char input[INPUTSIZE];
! 75: #define A_SIZE 1024
! 76: char a[A_SIZE];
! 77: int message = 0;
! 78: for (i=1; i<argc; i++) {
! 79: if (strcmp(argv[i],"-v") == 0) {
! 80: phc_verbose = 1;
! 81: }else if (strcmp(argv[i],"-g") == 0) {
! 82: phc_overwrite = 0;
! 83: }else if (strcmp(argv[i],"-i") == 0) {
! 84: ob = phc_call_phc(argv[i+1]);
! 85: n = phc_getoaSize(ob);
! 86: printf("[\n");
! 87: for (i=0; i<n; i++) {
! 88: phc_printObject(stdout,phc_getoa(ob,i));
! 89: if (i != n-1) printf(" ,\n"); else printf(" \n");
! 90: }
! 91: printf("]\n");
! 92: exit(0);
! 93: }
! 94: }
! 95: if (message) {
! 96: printf("Input example:\n 2 \n x**2 + y**2 - 1;\n x**2 + y**2 - 8*x - 3;\n");
! 97: printf("Note that input length is limited.\n");
! 98: }
! 99: while (1) {
! 100: if (message) printf("dim= ");
! 101: if (scanf("%d",&dim)<0) break;
! 102: sprintf(input,"%d\n",dim);
! 103: if (message) printf("Input %d equations please.\n",dim);
! 104: for (i=0; i<dim; i++) {
! 105: if (message) {printf("eq[%d] = ",i); fflush(stdout);}
! 106: do {
! 107: fgets(a,A_SIZE-1, stdin);
! 108: } while (strlen(a) == 0);
! 109: if (strlen(a) >= A_SIZE-3) {
! 110: fprintf(stderr,"Too big input for the input buffer a.\n"); exit(10);
! 111: }
! 112: if (strlen(input)+strlen(a) >= INPUTSIZE) {
! 113: fprintf(stderr,"Too big input for the input buffer input.\n"); exit(10);
! 114: }
! 115: sprintf(input+strlen(input),"%s\n",a);
! 116: }
! 117: ob = phc_call_phc(input);
! 118: if (message) {
! 119: printf("-----------------------------------------------------------\n");
! 120: }
! 121: n = phc_getoaSize(ob);
! 122: printf("[\n");
! 123: for (i=0; i<n; i++) {
! 124: phc_printObject(stdout,phc_getoa(ob,i));
! 125: if (i != n-1) printf(" ,\n"); else printf(" \n");
! 126: }
! 127: printf("]\n");
! 128: }
! 129: }
! 130:
! 131: int phc_scan_for_string(FILE *fp, char str[], int lenstr)
! 132: /*
! 133: ** Scans the file fp for a certain string str of length lenstr+1.
! 134: ** Reading stops when the string has been found, then the variable
! 135: ** on return equals 1, otherwise 0 is returned.
! 136: */
! 137: {
! 138: char buf[lenstr+1];
! 139: char ch;
! 140: int index,i,compare,npaths,dim,found;
! 141: index = -1;
! 142: found = 0;
! 143: while ((fscanf(fp,"%c",&ch)!=EOF) && found == 0)
! 144: {
! 145: if (index == -1 && ch == str[0])
! 146: {
! 147: index = 0;
! 148: buf[index] = ch;
! 149: }
! 150: else
! 151: {
! 152: if (index == lenstr)
! 153: {
! 154: compare = 0;
! 155: for (i=0; i<lenstr+1; i++)
! 156: {
! 157: if (buf[i]!=str[i])
! 158: {
! 159: compare = compare+1;
! 160: }
! 161: }
! 162: if (compare == 0)
! 163: {
! 164: found = 1;
! 165: }
! 166: index = -1;
! 167: }
! 168: else
! 169: if (index > -1 && index < lenstr)
! 170: {
! 171: index = index+1;
! 172: buf[index] = ch;
! 173: }
! 174: }
! 175: if (found == 1) break;
! 176: }
! 177: return found;
! 178: }
! 179: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim )
! 180: /*
! 181: ** Scans the file for the solutions, from a list of length npaths,
! 182: ** of complex vectors with dim entries.
! 183: ** The tolerance for the residual to a solution is set to 1.0E-12.
! 184: ** Returns solutions.
! 185: */
! 186: {
! 187: struct phc_object rob,sob;
! 188: char ch;
! 189: int fnd,i,j,nsols;
! 190: float res;
! 191: long double realpart;
! 192: long double imagpart;
! 193: long double realparts[npaths][dim];
! 194: long double imagparts[npaths][dim];
! 195: nsols = 0;
! 196: while (fscanf(fp,"%c",&ch)!=EOF)
! 197: {
! 198: fnd = phc_scan_for_string(fp,"start residual :",15);
! 199: if (fnd==1)
! 200: {
! 201: fscanf(fp,"%E",&res);
! 202: /* printf(" residual = "); printf("%E\n",res); */
! 203: if (res < 1.0E-12) nsols = nsols+1;
! 204: fnd = phc_scan_for_string(fp,"the solution for t :",19);
! 205: for (i=0;i<dim;i++)
! 206: {
! 207: fnd = phc_scan_for_string(fp,":",0);
! 208: fscanf(fp,"%LE",&realpart);
! 209: fscanf(fp,"%LE",&imagpart);
! 210: if (res < 1.0E-12)
! 211: {
! 212: realparts[nsols-1][i] = realpart;
! 213: imagparts[nsols-1][i] = imagpart;
! 214: }
! 215: }
! 216: }
! 217: }
! 218: if(phc_verbose) fprintf(stderr," number of solutions = %i\n",nsols);
! 219: rob = phc_newObjectArray(nsols);
! 220: for (i=0;i<nsols;i++)
! 221: {
! 222: /* fprintf(stderr,"Solution %i :\n",i+1); */
! 223: sob = phc_newObjectArray(dim);
! 224: for (j=0;j<dim;j++)
! 225: {
! 226: /*
! 227: printf("%20.14LE",realparts[i][j]); printf(" ");
! 228: printf("%20.14LE",imagparts[i][j]); printf("\n");
! 229: */
! 230: phc_putoa(sob,j,phc_complexTo(realparts[i][j],imagparts[i][j]));
! 231: }
! 232: phc_putoa(rob,i,sob);
! 233: }
! 234: return(rob);
! 235: }
! 236: struct phc_object phc_scan_output_of_phc(char *fname)
! 237: /*
! 238: ** Scans the file "output" of phc in two stages to get
! 239: ** 1) the number of paths and the dimension;
! 240: ** 2) the solutions, vectors with residuals < 1.0E-12.
! 241: */
! 242: {
! 243: FILE *otp;
! 244: char ch;
! 245: int fnd,npaths,dim,i,nsols;
! 246: otp = fopen(fname,"r");
! 247: if (phc_verbose) fprintf(stderr,"Scanning the %s of phc.\n",fname);
! 248: fnd = phc_scan_for_string(otp,"THE SOLUTIONS :",14);
! 249: fscanf(otp,"%i",&npaths);
! 250: if (phc_verbose) fprintf(stderr," number of paths traced = %i\n",npaths);
! 251: fscanf(otp,"%i",&dim);
! 252: if (phc_verbose) fprintf(stderr," dimension of solutions = %i\n",dim);
! 253: return(phc_scan_solutions(otp,npaths,dim));
! 254: }
! 255: struct phc_object phc_call_phc(char *sys) /* call phc, system as string */
! 256: {
! 257: FILE *inp;
! 258: char *f,*outf;
! 259: char cmd[1024];
! 260: char *w;
! 261: struct phc_object phc_NullObject ;
! 262: struct stat statbuf;
! 263:
! 264: phc_NullObject.tag = Snull;
! 265: f = phc_generateUniqueFileName("tmp.input");
! 266: if (phc_verbose) fprintf(stderr,"Creating file with name %s.\n",f);
! 267: outf = phc_generateUniqueFileName("tmp.output");
! 268: if (stat(outf,&statbuf) == 0) {
! 269: sprintf(cmd,"/bin/rm -f %s",outf);
! 270: system(cmd);
! 271: }
! 272: inp = fopen(f,"w");
! 273: fprintf(inp,sys);
! 274: fclose(inp);
! 275: if ((w = phc_which("phc")) != NULL) {
! 276: sprintf(cmd,"%s -b %s %s",w,f,outf);
! 277: }else{
! 278: sprintf(cmd,"phc -b %s %s",f,outf);
! 279: }
! 280: if (phc_verbose)fprintf(stderr,"Calling %s, black-box solver of phc.\n",cmd);
! 281: system(cmd);
! 282: if (stat(outf,&statbuf) < 0) {
! 283: fprintf(stderr,"execution error of phc.\n");
! 284: return(phc_NullObject);
! 285: }else{
! 286: if (phc_verbose) fprintf(stderr,"See the file %s for results.\n",outf);
! 287: return(phc_scan_output_of_phc(outf));
! 288: }
! 289: }
! 290:
! 291:
! 292: struct phc_object phc_newObjectArray(size)
! 293: int size;
! 294: {
! 295: struct phc_object rob;
! 296: struct phc_object *op;
! 297: if (size > 0) {
! 298: op = (struct phc_object *)sGC_malloc(size*sizeof(struct phc_object));
! 299: if (op == (struct phc_object *)NULL) {fprintf(stderr,"No memory\n");exit(1);}
! 300: }else{
! 301: op = (struct phc_object *)NULL;
! 302: }
! 303: rob.tag = Sarray;
! 304: rob.lc.ival = size;
! 305: rob.rc.op = op;
! 306: return(rob);
! 307: }
! 308:
! 309: void phc_printObject(FILE *fp,struct phc_object ob)
! 310: {
! 311: int n,i;
! 312: if (ob.tag == Snull) {
! 313: fprintf(fp,"null");
! 314: }else if (ob.tag == SstringObject) {
! 315: fprintf(fp,"%s",ob.lc.str);
! 316: }else if (ob.tag == Sarray) {
! 317: n = phc_getoaSize(ob);
! 318: fprintf(fp,"[");
! 319: for (i=0; i<n; i++) {
! 320: phc_printObject(fp,phc_getoa(ob,i));
! 321: if (i <n-1) fprintf(fp," , ");
! 322: }
! 323: fprintf(fp,"]");
! 324: }else if (ob.tag == SlongdoubleComplex) {
! 325: /* Try your favorite way to print complex numbers. */
! 326: /*fprintf(fp,"(%20.14LE)+I*(%20.14LE)",ob.lc.longdouble,ob.rc.longdouble);*/
! 327: fprintf(fp,"[%Lf , %Lf]",ob.lc.longdouble,ob.rc.longdouble);
! 328: }else{
! 329: fprintf(stderr,"Unknown phc_object tag %d",ob.tag);
! 330: }
! 331: }
! 332:
! 333:
! 334: char *phc_generateUniqueFileName(char *s)
! 335: {
! 336: char *t;
! 337: int i;
! 338: struct stat statbuf;
! 339: t = (char *)sGC_malloc(sizeof(char)*strlen(s)+4+2);
! 340: for (i=0; i<1000; i++) {
! 341: /* Give up if we failed for 1000 names. */
! 342: sprintf(t,"%s.%d",s,i);
! 343: if (phc_overwrite) return(t);
! 344: if (stat(t,&statbuf) < 0) {
! 345: return(t);
! 346: }
! 347: }
! 348: fprintf(stderr,"Could not generate a unique file name.");
! 349: return(NULL);
! 350: }
! 351:
! 352: char *phc_which(char *s) {
! 353: struct stat statbuf;
! 354: char *cmd,*a;
! 355: a = getenv("OpenXM_HOME");
! 356: if (a != NULL) {
! 357: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)+strlen(a)
! 358: +strlen("/usr/local/bin/")+1));
! 359: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
! 360: strcpy(cmd,a); strcat(cmd,"/bin/"); strcat(cmd,s);
! 361: if (stat(cmd,&statbuf) == 0) {
! 362: return(cmd);
! 363: }
! 364: }
! 365: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)
! 366: +strlen("/usr/local/bin/")+1));
! 367: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
! 368: strcpy(cmd,"/usr/local/bin/"); strcat(cmd,s);
! 369: if (stat(cmd,&statbuf) == 0) {
! 370: return(cmd);
! 371: }
! 372: strcpy(cmd,"/tmp/"); strcat(cmd,s);
! 373: if (stat(cmd,&statbuf) == 0) {
! 374: return(cmd);
! 375: }
! 376: return(NULL);
! 377: }
! 378:
! 379:
! 380: struct phc_object phc_complexTo(long double r, long double i)
! 381: {
! 382: struct phc_object rob;
! 383: rob.tag = SlongdoubleComplex;
! 384: rob.lc.longdouble = r;
! 385: rob.rc.longdouble = i;
! 386: return(rob);
! 387: }
! 388:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>