Annotation of OpenXM/src/phc/phc6.c, Revision 1.7
1.1 maekawa 1: /* phc6.c , yama:1999/sm1-prog/phc6.c */
1.7 ! takayama 2: /* $OpenXM: OpenXM/src/phc/phc6.c,v 1.6 2019/01/01 07:17:52 takayama Exp $ */
1.1 maekawa 3: /* This is a simple C interface to the black-box solver of phc.
4: ** Requirements:
5: ** 1) executable version of phc will be searched in the following order:
6: ** OpenXM_HOME/bin/, /usr/local/bin/phc, /tmp/phc, your command search path.
7: ** Here, PHC_LIBDIR is an environment variable.
8: ** 2) user of this program has write permissions to create
9: ** the files "tmp.input.xx" and "tmp.output.xx" in the directory where
10: ** this program is executed. xx are numbers.
11: */
12:
13: #include <stdio.h>
14: #include <sys/stat.h>
15: #include <unistd.h>
16: #include <stdlib.h>
1.6 takayama 17: #include <string.h>
18: #include "gc.h"
1.1 maekawa 19:
20: /* Definition of class identifiers. */
21: #define Snull 0
22: #define SstringObject 5
23: #define Sarray 6
24: #define SlongdoubleComplex 601
25:
26: /* Definition of Object */
27: union cell {
28: int ival;
29: char *str;
30: struct phc_object *op;
1.4 takayama 31: double longdouble;
1.1 maekawa 32: };
33: struct phc_object{
34: int tag; /* class identifier */
35: union cell lc; /* left cell */
36: union cell rc; /* right cell */
37: };
38:
39: /* Memory allocation function.
40: Use your favorite memory allocation function.
41: I recommend not to use malloc and to use gc4.14 for large applications. */
1.6 takayama 42: //#define sGC_malloc(n) malloc(n)
43: #define sGC_malloc(n) GC_malloc(n)
1.1 maekawa 44:
45: /********** macros to use Sarray **************/
46: /* put to Object Array */
47: #define phc_putoa(ob,i,cc) {\
1.4 takayama 48: if ((ob).tag != Sarray) {fprintf(stderr,"Warning: PUTOA is for an array of objects\n");} else \
49: {if ((0 <= (i)) && ((i) < (ob).lc.ival)) {\
50: (ob.rc.op)[i] = cc;\
51: }else{\
52: fprintf(stderr,"Warning: PUTOA, the size is %d.\n",(ob).lc.ival);\
53: }}}
1.1 maekawa 54: #define phc_getoa(ob,i) ((ob.rc.op)[i])
55: #define phc_getoaSize(ob) ((ob).lc.ival)
56:
57: /* prototypes */
58: struct phc_object phc_newObjectArray(int size);
59: void phc_printObject(FILE *fp,struct phc_object ob);
60: char *phc_generateUniqueFileName(char *s);
61: char *phc_which(char *s); /* search a path for the file s */
1.4 takayama 62: struct phc_object phc_complexTo(double r, double i);
1.1 maekawa 63:
64:
1.2 takayama 65: int phc_scan_for_string(FILE *fp, char str[]);
1.1 maekawa 66: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim );
1.6 takayama 67: struct phc_object phc_scan_output_of_phc(FILE *fp);
1.1 maekawa 68: struct phc_object phc_call_phc(char *sys);
69:
70: int phc_verbose = 0;
71: int phc_overwrite = 1; /* Always use tmp.input.0 and tmp.output.0
1.4 takayama 72: for work files. */
1.1 maekawa 73:
1.6 takayama 74: int main(int argc, char *argv[]) {
1.1 maekawa 75: struct phc_object ob;
1.6 takayama 76: struct phc_object ob2;
77: int n2;
1.1 maekawa 78: int n,i,dim;
79: #define INPUTSIZE 4096
80: char input[INPUTSIZE];
1.4 takayama 81: char fname[INPUTSIZE];
1.1 maekawa 82: #define A_SIZE 1024
83: char a[A_SIZE];
84: int message = 0;
1.6 takayama 85: FILE *fp;
1.1 maekawa 86: for (i=1; i<argc; i++) {
87: if (strcmp(argv[i],"-v") == 0) {
1.2 takayama 88: phc_verbose = 1; message=1;
1.1 maekawa 89: }else if (strcmp(argv[i],"-g") == 0) {
90: phc_overwrite = 0;
1.4 takayama 91: }else if (strcmp(argv[i],"-file") == 0) {
92: /* For debugging. */
93: i++;
94: strncpy(fname,argv[i],INPUTSIZE-1);
1.6 takayama 95: fp = fopen(fname,"r");
96: if (fp == NULL) {
97: fprintf(stderr,"File %s is not found.\n",fname);
98: return(-1);
99: }
100: ob = phc_scan_output_of_phc(fp); /* stable mixed vol, out of torus */
1.4 takayama 101: n = phc_getoaSize(ob);
1.6 takayama 102: ob2 = phc_scan_output_of_phc(fp); /* in torus */
103: n2 = phc_getoaSize(ob2);
104: if (phc_verbose) fprintf(stderr,"n=%d, n2=%d\n",n,n2);
1.4 takayama 105: printf("[\n");
106: for (i=0; i<n; i++) {
107: phc_printObject(stdout,phc_getoa(ob,i));
1.6 takayama 108: if ((i != n-1) || (n2 > 0)) printf(" ,\n"); else printf(" \n");
109: }
110: for (i=0; i<n2; i++) {
111: phc_printObject(stdout,phc_getoa(ob2,i));
112: if (i != n2-1) printf(" ,\n"); else printf(" \n");
1.4 takayama 113: }
114: printf("]\n");
1.6 takayama 115: return(0);
1.1 maekawa 116: }else if (strcmp(argv[i],"-i") == 0) {
117: ob = phc_call_phc(argv[i+1]);
118: n = phc_getoaSize(ob);
119: printf("[\n");
120: for (i=0; i<n; i++) {
1.4 takayama 121: phc_printObject(stdout,phc_getoa(ob,i));
122: if (i != n-1) printf(" ,\n"); else printf(" \n");
1.1 maekawa 123: }
124: printf("]\n");
125: exit(0);
126: }
127: }
128: if (message) {
129: printf("Input example:\n 2 \n x**2 + y**2 - 1;\n x**2 + y**2 - 8*x - 3;\n");
130: printf("Note that input length is limited.\n");
131: }
132: while (1) {
133: if (message) printf("dim= ");
1.4 takayama 134: if (fgets(input,INPUTSIZE,stdin) <= 0) break;
1.2 takayama 135: sscanf(input,"%d",&dim);
1.1 maekawa 136: sprintf(input,"%d\n",dim);
137: if (message) printf("Input %d equations please.\n",dim);
138: for (i=0; i<dim; i++) {
1.4 takayama 139: if (message) {printf("eq[%d] = ",i); fflush(stdout);}
140: do {
141: fgets(a,A_SIZE-1, stdin);
142: } while (strlen(a) == 0);
143: if (strlen(a) >= A_SIZE-3) {
144: fprintf(stderr,"Too big input for the input buffer a.\n"); exit(10);
145: }
146: if (strlen(input)+strlen(a) >= INPUTSIZE) {
147: fprintf(stderr,"Too big input for the input buffer input.\n"); exit(10);
148: }
149: sprintf(input+strlen(input),"%s\n",a);
1.1 maekawa 150: }
151: ob = phc_call_phc(input);
152: if (message) {
153: printf("-----------------------------------------------------------\n");
154: }
155: n = phc_getoaSize(ob);
156: printf("[\n");
157: for (i=0; i<n; i++) {
158: phc_printObject(stdout,phc_getoa(ob,i));
159: if (i != n-1) printf(" ,\n"); else printf(" \n");
160: }
161: printf("]\n");
162: }
1.6 takayama 163: return(0);
1.1 maekawa 164: }
165:
1.4 takayama 166: int phc_scan_for_string(FILE *fp, char s[])
167: /*
168: ** Scans the file fp for a certain string str of length lenstr+1.
169: ** Reading stops when the string has been found, then the variable
170: ** on return equals 1, otherwise 0 is returned.
171: */
172: #define BUF_SIZE 1024
173: {
174: char buf[BUF_SIZE];
175: int pt,n,c,i;
176: pt = 0;
177: n=strlen(s);
178: if (n > BUF_SIZE-2) {
179: fprintf(stderr,"Too long string for scan_for_string.\n");
180: exit(1);
181: }
182: for (i=0; i<n; i++) {
183: buf[i] = fgetc(fp); buf[i+1]='\0';
184: if (buf[i] == EOF) return 0;
185: }
1.6 takayama 186: if (strcmp(s,buf) == 0) return 1;
1.4 takayama 187: while ((c=fgetc(fp)) != EOF) {
188: for (i=0; i<n; i++) {
189: buf[i] = buf[i+1];
190: }
191: buf[n-1]=c; buf[n] = '\0';
192: if (strcmp(s,buf) == 0) return 1;
193: }
194: return 0;
195: }
196:
197: int phc_scan_for_string_old(FILE *fp, char str[])
1.1 maekawa 198: /*
199: ** Scans the file fp for a certain string str of length lenstr+1.
200: ** Reading stops when the string has been found, then the variable
201: ** on return equals 1, otherwise 0 is returned.
202: */
203: {
1.2 takayama 204: #define BUF_SIZE 1024
205: char buf[BUF_SIZE];
1.1 maekawa 206: char ch;
207: int index,i,compare,npaths,dim,found;
1.2 takayama 208: int lenstr;
209: lenstr = strlen(str);
210: if (lenstr >= BUF_SIZE-1) {
1.4 takayama 211: fprintf(stderr,"Too long string in phc_scan_for_string\n");
212: exit(-1);
1.2 takayama 213: }
1.1 maekawa 214: index = -1;
215: found = 0;
1.2 takayama 216: while (((ch = fgetc(fp))!=EOF) && found == 0)
1.4 takayama 217: {
218: if (index == -1 && ch == str[0])
219: {
220: index = 0;
221: buf[index] = ch;
222: }
223: else
224: {
225: if (index == lenstr)
226: {
227: compare = 0;
228: for (i=0; str[i] != '\0'; i++)
229: {
230: if (buf[i]!=str[i])
231: {
232: compare = compare+1;
233: }
234: }
235: if (compare == 0)
236: {
237: found = 1;
238: }
239: index = -1;
240: }
241: else
242: if (index > -1 && index < lenstr)
243: {
244: index = index+1;
245: buf[index] = ch;
246: }
247: }
248: if (found == 1) break;
249: }
1.1 maekawa 250: return found;
251: }
252: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim )
253: /*
254: ** Scans the file for the solutions, from a list of length npaths,
255: ** of complex vectors with dim entries.
256: ** The tolerance for the residual to a solution is set to 1.0E-12.
257: ** Returns solutions.
258: */
1.4 takayama 259: #define BUFSIZE 1024
1.1 maekawa 260: {
261: struct phc_object rob,sob;
262: char ch;
263: int fnd,i,j,nsols;
1.2 takayama 264: double res;
1.4 takayama 265: double realpart;
266: double imagpart;
267: /* double realparts[npaths][dim];
268: double imagparts[npaths][dim]; */
269: double *realparts;
270: double *imagparts;
1.6 takayama 271: double *res_list;
1.4 takayama 272: char buf[BUFSIZE];
1.1 maekawa 273: nsols = 0;
1.5 takayama 274: realparts = (double *)sGC_malloc(sizeof(double)*(npaths*dim+1));
275: imagparts = (double *)sGC_malloc(sizeof(double)*(npaths*dim+1));
1.6 takayama 276: res_list = (double *)sGC_malloc(sizeof(double)*npaths+1);
1.2 takayama 277: while ((ch = fgetc(fp)) != EOF)
1.1 maekawa 278: {
1.2 takayama 279: fnd = phc_scan_for_string(fp,"start residual :");
1.1 maekawa 280: if (fnd==1)
1.4 takayama 281: {
1.6 takayama 282: fgets(buf,BUFSIZE-1,fp);
1.7 ! takayama 283: sscanf(buf,"%lg",&res);
1.6 takayama 284: if (phc_verbose) {
285: fprintf(stderr,"res in string: %s\n",buf);
286: fprintf(stderr," residual = "); fprintf(stderr,"%le\n",res);
287: }
288: nsols++;
289: /*
1.4 takayama 290: if (res < 1.0E-12) {
291: nsols = nsols+1;
292: if (nsols > npaths) {
293: fprintf(stderr,"Something is wrong in phc_scan_solutions\n");
294: fprintf(stderr,"npaths=%d, nsols=%d \n",npaths,nsols);
295: exit(-1);
296: }
297: }
1.6 takayama 298: */
1.4 takayama 299: fnd = phc_scan_for_string(fp,"the solution for t :");
1.6 takayama 300: res_list[nsols-1] = res;
1.4 takayama 301: for (i=0;i<dim;i++)
302: {
303: fnd = phc_scan_for_string(fp,":");
1.6 takayama 304: fgets(buf,BUFSIZE-1,fp);
1.7 ! takayama 305: sscanf(buf,"%lf %lf",&realpart,&imagpart);
1.6 takayama 306: if (phc_verbose) fprintf(stderr,"sols in string: %s\n",buf);
1.4 takayama 307: /*fscanf(fp,"%le",&realpart);
1.6 takayama 308: fscanf(fp,"%le",&imagpart);*/
309: realparts[(nsols-1)*dim+i] = realpart;
310: imagparts[(nsols-1)*dim+i] = imagpart;
1.4 takayama 311: }
1.6 takayama 312: if (nsols >= npaths) break;
1.4 takayama 313: }
1.1 maekawa 314: }
1.6 takayama 315: if (nsols < npaths) {
316: fprintf(stderr,"Warning: nsols < npaths. Check tmp.output.*\n");
317: }
1.1 maekawa 318: if(phc_verbose) fprintf(stderr," number of solutions = %i\n",nsols);
319: rob = phc_newObjectArray(nsols);
320: for (i=0;i<nsols;i++)
321: {
1.5 takayama 322: if (phc_verbose) fprintf(stderr,"Solution %i :\n",i+1);
1.1 maekawa 323: sob = phc_newObjectArray(dim);
324: for (j=0;j<dim;j++)
1.4 takayama 325: {
326: /*
327: printf("%20.14le",realparts[i][j]); printf(" ");
328: printf("%20.14le",imagparts[i][j]); printf("\n");
329: */
1.5 takayama 330: if (phc_verbose) {
331: printf("%20.14le",realparts[i*dim+j]); printf(" ");
332: printf("%20.14le",imagparts[i*dim+j]); printf("\n");
333: }
1.4 takayama 334:
335: /*phc_putoa(sob,j,phc_complexTo(realparts[i][j],imagparts[i][j]));*/
1.5 takayama 336: phc_putoa(sob,j,phc_complexTo(realparts[i*dim+j],
337: imagparts[i*dim+j]));
1.4 takayama 338: }
1.1 maekawa 339: phc_putoa(rob,i,sob);
340: }
1.6 takayama 341: /* Todo, res_list[] is not returned. */
1.1 maekawa 342: return(rob);
343: }
1.6 takayama 344: struct phc_object phc_scan_output_of_phc(FILE *fp)
1.1 maekawa 345: /*
346: ** Scans the file "output" of phc in two stages to get
347: ** 1) the number of paths and the dimension;
348: ** 2) the solutions, vectors with residuals < 1.0E-12.
349: */
350: {
351: FILE *otp;
352: char ch;
353: int fnd,npaths,dim,i,nsols;
1.6 takayama 354: struct phc_object rob;
355: otp = fp;
1.4 takayama 356: if (otp == NULL) {
1.6 takayama 357: fprintf(stderr,"Error: the file is not found.\n");
1.4 takayama 358: exit(0);
359: }
1.2 takayama 360: fnd = phc_scan_for_string(otp,"THE SOLUTIONS :");
1.6 takayama 361: if (fnd < 1) {
362: if (phc_verbose) fprintf(stderr,"No more THE SOLUTIONS\n");
363: rob = phc_newObjectArray(0);
364: return(rob);
365: }
1.1 maekawa 366: fscanf(otp,"%i",&npaths);
367: if (phc_verbose) fprintf(stderr," number of paths traced = %i\n",npaths);
368: fscanf(otp,"%i",&dim);
369: if (phc_verbose) fprintf(stderr," dimension of solutions = %i\n",dim);
370: return(phc_scan_solutions(otp,npaths,dim));
371: }
372: struct phc_object phc_call_phc(char *sys) /* call phc, system as string */
373: {
374: FILE *inp;
375: char *f,*outf;
376: char cmd[1024];
377: char *w;
378: struct phc_object phc_NullObject ;
379: struct stat statbuf;
1.6 takayama 380: FILE *fp;
381: struct phc_object ob;
382: struct phc_object ob2;
383: struct phc_object ob_all;
384: int i,n,n2;
1.1 maekawa 385:
386: phc_NullObject.tag = Snull;
387: f = phc_generateUniqueFileName("tmp.input");
388: if (phc_verbose) fprintf(stderr,"Creating file with name %s.\n",f);
389: outf = phc_generateUniqueFileName("tmp.output");
390: if (stat(outf,&statbuf) == 0) {
391: sprintf(cmd,"/bin/rm -f %s",outf);
392: system(cmd);
393: }
394: inp = fopen(f,"w");
1.6 takayama 395: fprintf(inp,"%s",sys);
1.1 maekawa 396: fclose(inp);
397: if ((w = phc_which("phc")) != NULL) {
398: sprintf(cmd,"%s -b %s %s",w,f,outf);
399: }else{
400: sprintf(cmd,"phc -b %s %s",f,outf);
401: }
402: if (phc_verbose)fprintf(stderr,"Calling %s, black-box solver of phc.\n",cmd);
403: system(cmd);
404: if (stat(outf,&statbuf) < 0) {
405: fprintf(stderr,"execution error of phc.\n");
406: return(phc_NullObject);
407: }else{
408: if (phc_verbose) fprintf(stderr,"See the file %s for results.\n",outf);
1.6 takayama 409: fp = fopen(outf,"r");
410: ob=phc_scan_output_of_phc(fp);
411: n = phc_getoaSize(ob);
412: ob2=phc_scan_output_of_phc(fp);
413: n2 = phc_getoaSize(ob2);
414: ob_all=phc_newObjectArray(n+n2);
415: if (phc_verbose) fprintf(stderr,"n=%d, n2=%d\n",n,n2);
416: for (i=0; i<n; i++) phc_putoa(ob_all,i,phc_getoa(ob,i));
417: for (i=0; i<n2; i++) phc_putoa(ob_all,n+i,phc_getoa(ob2,i));
418: return(ob_all);
1.1 maekawa 419: }
420: }
421:
422:
423: struct phc_object phc_newObjectArray(size)
1.4 takayama 424: int size;
1.1 maekawa 425: {
426: struct phc_object rob;
427: struct phc_object *op;
428: if (size > 0) {
429: op = (struct phc_object *)sGC_malloc(size*sizeof(struct phc_object));
430: if (op == (struct phc_object *)NULL) {fprintf(stderr,"No memory\n");exit(1);}
431: }else{
432: op = (struct phc_object *)NULL;
433: }
434: rob.tag = Sarray;
435: rob.lc.ival = size;
436: rob.rc.op = op;
437: return(rob);
438: }
439:
440: void phc_printObject(FILE *fp,struct phc_object ob)
441: {
442: int n,i;
443: if (ob.tag == Snull) {
444: fprintf(fp,"null");
445: }else if (ob.tag == SstringObject) {
446: fprintf(fp,"%s",ob.lc.str);
447: }else if (ob.tag == Sarray) {
448: n = phc_getoaSize(ob);
449: fprintf(fp,"[");
450: for (i=0; i<n; i++) {
451: phc_printObject(fp,phc_getoa(ob,i));
452: if (i <n-1) fprintf(fp," , ");
453: }
454: fprintf(fp,"]");
455: }else if (ob.tag == SlongdoubleComplex) {
456: /* Try your favorite way to print complex numbers. */
1.4 takayama 457: /*fprintf(fp,"(%20.14le)+I*(%20.14le)",ob.lc.longdouble,ob.rc.longdouble);*/
1.7 ! takayama 458: fprintf(fp,"[%.16lg , %.16lg]",ob.lc.longdouble,ob.rc.longdouble);
1.1 maekawa 459: }else{
460: fprintf(stderr,"Unknown phc_object tag %d",ob.tag);
461: }
462: }
463:
464:
465: char *phc_generateUniqueFileName(char *s)
466: {
467: char *t;
468: int i;
469: struct stat statbuf;
470: t = (char *)sGC_malloc(sizeof(char)*strlen(s)+4+2);
471: for (i=0; i<1000; i++) {
472: /* Give up if we failed for 1000 names. */
473: sprintf(t,"%s.%d",s,i);
474: if (phc_overwrite) return(t);
475: if (stat(t,&statbuf) < 0) {
476: return(t);
477: }
478: }
479: fprintf(stderr,"Could not generate a unique file name.");
480: return(NULL);
481: }
482:
483: char *phc_which(char *s) {
484: struct stat statbuf;
485: char *cmd,*a;
486: a = getenv("OpenXM_HOME");
487: if (a != NULL) {
488: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)+strlen(a)
1.4 takayama 489: +strlen("/usr/local/bin/")+1));
1.1 maekawa 490: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
491: strcpy(cmd,a); strcat(cmd,"/bin/"); strcat(cmd,s);
492: if (stat(cmd,&statbuf) == 0) {
493: return(cmd);
494: }
495: }
496: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)
1.4 takayama 497: +strlen("/usr/local/bin/")+1));
1.1 maekawa 498: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
499: strcpy(cmd,"/usr/local/bin/"); strcat(cmd,s);
500: if (stat(cmd,&statbuf) == 0) {
501: return(cmd);
502: }
503: strcpy(cmd,"/tmp/"); strcat(cmd,s);
504: if (stat(cmd,&statbuf) == 0) {
505: return(cmd);
506: }
507: return(NULL);
508: }
509:
510:
1.4 takayama 511: struct phc_object phc_complexTo(double r, double i)
1.1 maekawa 512: {
513: struct phc_object rob;
514: rob.tag = SlongdoubleComplex;
515: rob.lc.longdouble = r;
516: rob.rc.longdouble = i;
517: return(rob);
518: }
519:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to phc6.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / phc