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