Annotation of OpenXM/src/ox_gsl/ox_gsl.c, Revision 1.15
1.15 ! takayama 1: /* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.14 2018/06/07 11:49:51 takayama Exp $
1.1 takayama 2: */
3:
4: #include <stdio.h>
5: #include <stdlib.h>
1.2 takayama 6: #include <setjmp.h>
7: #include <string.h>
1.3 takayama 8: #include <unistd.h>
1.10 ohara 9: #include <signal.h>
1.3 takayama 10: #include <math.h>
11: #include "ox_gsl.h"
12: #include "call_gsl.h" // need only when you bind call_gsl functions.
1.1 takayama 13:
14: OXFILE *fd_rw;
15:
16: #define INIT_S_SIZE 2048
17: #define EXT_S_SIZE 2048
18:
19: static int stack_size = 0;
20: static int stack_pointer = 0;
21: static cmo **stack = NULL;
22:
23: int Debug=1;
24:
1.2 takayama 25: void show_stack_top() {
1.1 takayama 26: cmo *data;
27: if (stack_pointer > 0) {
28: data=stack[stack_pointer-1];
1.2 takayama 29: print_cmo(data); printf("\n");
1.1 takayama 30: }else {
31: printf("The stack is empty.\n");
32: }
33: }
34:
1.3 takayama 35: void *gc_realloc(void *p,size_t osize,size_t nsize)
36: { return (void *)GC_realloc(p,nsize);}
37:
38: void gc_free(void *p,size_t size)
39: { GC_free(p);}
40:
1.1 takayama 41: void init_gc()
1.3 takayama 42: { GC_INIT();
43: mp_set_memory_functions(GC_malloc,gc_realloc,gc_free);
1.5 takayama 44: init_dic(); // initialize ox_eval.c
1.1 takayama 45: }
46:
1.2 takayama 47: void initialize_stack()
1.1 takayama 48: {
49: stack_pointer = 0;
50: stack_size = INIT_S_SIZE;
1.3 takayama 51: stack = GC_malloc(stack_size*sizeof(cmo*));
1.1 takayama 52: }
53:
1.2 takayama 54: static void extend_stack()
1.1 takayama 55: {
56: int size2 = stack_size + EXT_S_SIZE;
57: cmo **stack2 = malloc(size2*sizeof(cmo*));
58: memcpy(stack2, stack, stack_size*sizeof(cmo *));
59: free(stack);
60: stack = stack2;
61: stack_size = size2;
62: }
63:
1.2 takayama 64: void push(cmo* m)
1.1 takayama 65: {
66: stack[stack_pointer] = m;
67: stack_pointer++;
68: if (stack_pointer >= stack_size) {
69: extend_stack();
70: }
71: }
72:
73: cmo* pop()
74: {
75: if (stack_pointer > 0) {
76: stack_pointer--;
77: return stack[stack_pointer];
78: }
79: return new_cmo_null();
80: }
81:
82: void pops(int n)
83: {
84: stack_pointer -= n;
85: if (stack_pointer < 0) {
86: stack_pointer = 0;
87: }
88: }
89:
90: #define OX_GSL_VERSION 2018032901
91: #define ID_STRING "2018/03/29 13:56:00"
92:
93: int sm_mathcap()
94: {
1.2 takayama 95: int available_cmo[]={
96: CMO_NULL,
97: CMO_INT32,
98: // CMO_DATUM,
99: CMO_STRING,
100: CMO_MATHCAP,
101: CMO_LIST,
102: // CMO_MONOMIAL32,
103: CMO_ZZ,
104: // CMO_QQ,
105: CMO_BIGFLOAT32,
106: CMO_COMPLEX,
107: CMO_IEEE_DOUBLE_FLOAT,
108: CMO_ZERO,
109: // CMO_DMS_GENERIC,
110: // CMO_RING_BY_NAME,
111: // CMO_INDETERMINATE,
112: // CMO_DISTRIBUTED_POLYNOMIAL,
113: // CMO_RECURSIVE_POLYNOMIAL,
114: // CMO_POLYNOMIAL_IN_ONE_VARIABLE,
1.5 takayama 115: CMO_TREE,
1.2 takayama 116: CMO_ERROR2,
117: 0};
118: int available_sm_command[]={
119: SM_popCMO,
120: SM_popString,
121: SM_mathcap,
122: SM_pops,
123: // SM_executeStringByLocalParser,
124: SM_executeFunction,
125: SM_setMathCap,
126: SM_shutdown,
127: SM_control_kill,
128: SM_control_reset_connection,
129: SM_control_spawn_server,
130: SM_control_terminate_server,
131: 0};
132: mathcap_init(OX_GSL_VERSION, ID_STRING, "ox_gsl", available_cmo,available_sm_command);
133: push((cmo *)oxf_cmo_mathcap(fd_rw));
1.1 takayama 134: return 0;
135: }
136:
137: int sm_popCMO()
138: {
139: cmo* m = pop();
140:
141: if (m != NULL) {
142: send_ox_cmo(fd_rw, m);
143: return 0;
144: }
145: return SM_popCMO;
146: }
147:
1.3 takayama 148: cmo *make_error2(const char *reason,const char *fname,int line,int code)
1.1 takayama 149: {
1.3 takayama 150: // gsl_error_handler_t void handler(const char *reason,const char *file,int line, int gsl_errno)
151: cmo *ms;
152: cmo *err;
153: cmo *m;
154: cmo **argv;
155: int n;
156: char *s;
157: n = 5;
158: argv = (cmo **) GC_malloc(sizeof(cmo *)*n);
159: ms = (cmo *)new_cmo_string("Error"); argv[0] = ms;
1.4 takayama 160: if (reason != NULL) {s = (char *)GC_malloc(strlen(reason)+1); strcpy(s,reason);
161: }else strcpy(s,"");
1.3 takayama 162: ms = (cmo *) new_cmo_string(s); argv[1] = ms;
1.4 takayama 163: if (fname != NULL) {s = (char *)GC_malloc(strlen(fname)+1); strcpy(s,fname);
164: }else strcpy(s,"");
1.3 takayama 165: ms = (cmo *) new_cmo_string(s); argv[2] = ms;
166: err = (cmo *)new_cmo_int32(line); argv[3] = err;
167: err = (cmo *)new_cmo_int32(code); argv[4] = err;
168:
169: m = (cmo *)new_cmo_list_array((void *)argv,n);
170: return (m);
1.1 takayama 171: }
172:
173: int get_i()
174: {
175: cmo *c = pop();
176: if (c->tag == CMO_INT32) {
177: return ((cmo_int32 *)c)->i;
178: }else if (c->tag == CMO_ZZ) {
179: return mpz_get_si(((cmo_zz *)c)->mpz);
1.2 takayama 180: }else if (c->tag == CMO_NULL) {
181: return(0);
182: }else if (c->tag == CMO_ZERO) {
183: return(0);
1.1 takayama 184: }
1.4 takayama 185: myhandler("get_i: not an integer",NULL,0,-1);
1.1 takayama 186: return 0;
187: }
188:
1.2 takayama 189: void get_xy(int *x, int *y)
1.1 takayama 190: {
191: pop();
192: *x = get_i();
193: *y = get_i();
194: }
195:
1.2 takayama 196: void my_add_int32()
1.1 takayama 197: {
198: int x, y;
199: get_xy(&x, &y);
1.2 takayama 200: push((cmo *)new_cmo_int32(x+y));
1.1 takayama 201: }
202:
1.2 takayama 203: double get_double()
204: {
1.3 takayama 205: #define mympz(c) (((cmo_zz *)c)->mpz)
1.2 takayama 206: cmo *c = pop();
207: if (c->tag == CMO_INT32) {
208: return( (double) (((cmo_int32 *)c)->i) );
209: }else if (c->tag == CMO_IEEE_DOUBLE_FLOAT) {
1.3 takayama 210: return (((cmo_double *)c)->d); // see ox_toolkit.h
1.2 takayama 211: }else if (c->tag == CMO_ZZ) {
1.3 takayama 212: if ((mpz_cmp_si(mympz(c),(long int) 0x7fffffff)>0) ||
213: (mpz_cmp_si(mympz(c),(long int) -0x7fffffff)<0)) {
1.4 takayama 214: myhandler("get_double: out of int32",NULL,0,-1);
1.3 takayama 215: return(NAN);
216: }
217: return( (double) mpz_get_si(((cmo_zz *)c)->mpz));
1.2 takayama 218: }else if (c->tag == CMO_NULL) {
219: return(0);
220: }else if (c->tag == CMO_ZERO) {
221: return(0);
222: }
1.4 takayama 223: myhandler("get_double: not a double",NULL,0,-1);
1.3 takayama 224: return(NAN);
1.2 takayama 225: }
1.14 takayama 226: /* get_double() will be obsolted and will be replaced by cmo2double(c) */
227: double cmo2double(cmo *c)
228: {
229: #define mympz(c) (((cmo_zz *)c)->mpz)
230: if (c == NULL) c = pop();
231: if (c->tag == CMO_INT32) {
232: return( (double) (((cmo_int32 *)c)->i) );
233: }else if (c->tag == CMO_IEEE_DOUBLE_FLOAT) {
234: return (((cmo_double *)c)->d); // see ox_toolkit.h
235: }else if (c->tag == CMO_ZZ) {
236: if ((mpz_cmp_si(mympz(c),(long int) 0x7fffffff)>0) ||
237: (mpz_cmp_si(mympz(c),(long int) -0x7fffffff)<0)) {
238: myhandler("get_double: out of int32",NULL,0,-1);
239: return(NAN);
240: }
241: return( (double) mpz_get_si(((cmo_zz *)c)->mpz));
242: }else if (c->tag == CMO_NULL) {
243: return(0);
244: }else if (c->tag == CMO_ZERO) {
245: return(0);
246: }
247: myhandler("cmo2double: not a double",NULL,0,-1);
248: return(NAN);
249: }
1.2 takayama 250:
251: void my_add_double() {
252: double x,y;
253: pop();
254: y = get_double();
255: x = get_double();
256: push((cmo *)new_cmo_double(x+y));
257: }
258:
259: double *get_double_list(int *length) {
260: cmo *c;
261: cmo *entry;
262: cell *cellp;
263: double *d;
264: int n,i;
265: c = pop();
266: if (c->tag != CMO_LIST) {
1.3 takayama 267: // make_error2("get_double_list",NULL,0,-1);
1.2 takayama 268: *length=-1; return(0);
269: }
270: n = *length = list_length((cmo_list *)c);
271: d = (double *) GC_malloc(sizeof(double)*(*length+1));
272: cellp = list_first((cmo_list *)c);
273: entry = cellp->cmo;
274: for (i=0; i<n; i++) {
275: if (Debug) {
276: printf("entry[%d]=",i); print_cmo(entry); printf("\n");
277: }
278: if (entry->tag == CMO_INT32) {
279: d[i]=( (double) (((cmo_int32 *)entry)->i) );
280: }else if (entry->tag == CMO_IEEE_DOUBLE_FLOAT) {
281: d[i]=((cmo_double *)entry)->d;
282: }else if (entry->tag == CMO_ZZ) {
283: d[i]=( (double) mpz_get_si(((cmo_zz *)entry)->mpz));
284: }else if (entry->tag == CMO_NULL) {
285: d[i]= 0;
286: }else {
287: fprintf(stderr,"entries of the list should be int32 or zz or double\n");
288: *length = -1;
1.3 takayama 289: myhandler("get_double_list",NULL,0,-1);
1.2 takayama 290: return(NULL);
291: }
292: cellp = list_next(cellp);
293: entry = cellp->cmo;
294: }
295: return(d);
296: }
1.14 takayama 297: /* get_double_list will be obsolted and will be replaced by cmo2double_list() */
298: double *cmo2double_list(int *length,cmo *c) {
299: cmo *entry;
300: cell *cellp;
301: double *d;
302: int n,i;
303: if (c == NULL) c = pop();
304: if (c->tag != CMO_LIST) {
305: // make_error2("get_double_list",NULL,0,-1);
306: *length=-1; return(0);
307: }
308: n = *length = list_length((cmo_list *)c);
309: d = (double *) GC_malloc(sizeof(double)*(*length+1));
310: cellp = list_first((cmo_list *)c);
311: entry = cellp->cmo;
312: for (i=0; i<n; i++) {
313: if (Debug) {
314: printf("entry[%d]=",i); print_cmo(entry); printf("\n");
315: }
316: if (entry->tag == CMO_INT32) {
317: d[i]=( (double) (((cmo_int32 *)entry)->i) );
318: }else if (entry->tag == CMO_IEEE_DOUBLE_FLOAT) {
319: d[i]=((cmo_double *)entry)->d;
320: }else if (entry->tag == CMO_ZZ) {
321: d[i]=( (double) mpz_get_si(((cmo_zz *)entry)->mpz));
322: }else if (entry->tag == CMO_NULL) {
323: d[i]= 0;
324: }else {
325: fprintf(stderr,"entries of the list should be int32 or zz or double\n");
326: *length = -1;
327: myhandler("get_double_list",NULL,0,-1);
328: return(NULL);
329: }
330: cellp = list_next(cellp);
331: entry = cellp->cmo;
332: }
333: return(d);
334: }
1.2 takayama 335: void show_double_list() {
336: int n;
337: double *d;
338: int i;
339: pop(); // pop argument number;
340: d = get_double_list(&n);
1.3 takayama 341: if (n < 0) fprintf(stderr,"Error in the double list\n");
1.2 takayama 342: printf("show_double_list: length=%d\n",n);
343: for (i=0; i<n; i++) {
344: printf("%lg, ",d[i]);
345: }
346: printf("\n");
347: }
348:
349: char *get_string() {
350: cmo *c;
351: c = pop();
352: if (c->tag == CMO_STRING) {
353: return (((cmo_string *)c)->s);
354: }
1.3 takayama 355: // make_error2(-1);
1.2 takayama 356: return(NULL);
357: }
1.1 takayama 358:
1.7 takayama 359: void test_ox_eval() {
1.5 takayama 360: cmo *c;
361: double d=0;
362: pop();
1.7 takayama 363: c=pop();
1.5 takayama 364: if (Debug) {
1.7 takayama 365: ox_printf("cmo *c="); print_cmo(c); ox_printf("\n");
1.5 takayama 366: }
1.6 ohara 367: init_dic();
1.5 takayama 368: register_entry("x",1.25);
1.7 takayama 369: if (eval_cmo(c,&d) == 0) myhandler("eval_cmo failed",NULL,0,-1);
1.5 takayama 370: push((cmo *)new_cmo_double(d));
371: }
372:
1.1 takayama 373: int sm_executeFunction()
374: {
375: cmo_string *func = (cmo_string *)pop();
376: if (func->tag != CMO_STRING) {
1.3 takayama 377: push(make_error2("sm_executeFunction, not CMO_STRING",NULL,0,-1));
1.1 takayama 378: return -1;
379: }
1.9 takayama 380: init_dic();
1.2 takayama 381: // Test functions
1.1 takayama 382: if (strcmp(func->s, "add_int32") == 0) {
383: my_add_int32();
1.2 takayama 384: }else if (strcmp(func->s,"add_double")==0) {
385: my_add_double();
386: }else if (strcmp(func->s,"show_double_list")==0) {
387: show_double_list();
1.3 takayama 388: }else if (strcmp(func->s,"restart")==0) {
389: pop(); restart();
1.5 takayama 390: }else if (strcmp(func->s,"test_ox_eval")==0) {
391: test_ox_eval();
1.2 takayama 392: // The following functions are defined in call_gsl.c
393: }else if (strcmp(func->s,"gsl_sf_lngamma_complex_e")==0) {
394: call_gsl_sf_lngamma_complex_e();
1.8 takayama 395: }else if (strcmp(func->s,"gsl_integration_qags")==0) {
396: call_gsl_integration_qags();
1.11 takayama 397: }else if (strcmp(func->s,"gsl_monte_plain_integrate")==0) {
1.12 takayama 398: call_gsl_monte_plain_miser_vegas_integrate(0);
399: }else if (strcmp(func->s,"gsl_monte_miser_integrate")==0) {
400: call_gsl_monte_plain_miser_vegas_integrate(1);
401: }else if (strcmp(func->s,"gsl_monte_vegas_integrate")==0) {
402: call_gsl_monte_plain_miser_vegas_integrate(2);
1.15 ! takayama 403: }else if (strcmp(func->s,"gsl_odeiv_step_rk4")==0) {
! 404: call_gsl_odeiv_step("rk4");
1.1 takayama 405: }else {
1.3 takayama 406: push(make_error2("sm_executeFunction, unknown function",NULL,0,-1));
1.1 takayama 407: return -1;
1.2 takayama 408: }
409: return(0);
1.1 takayama 410: }
411:
412:
413: int receive_and_execute_sm_command()
414: {
415: int code = receive_int32(fd_rw);
416: switch(code) {
417: case SM_popCMO:
418: sm_popCMO();
419: break;
420: case SM_executeFunction:
421: sm_executeFunction();
422: break;
423: case SM_mathcap:
424: sm_mathcap();
425: break;
426: case SM_setMathCap:
427: pop();
428: break;
429: default:
430: ;
431: }
1.2 takayama 432: return(0);
1.1 takayama 433: }
434:
435: int receive()
436: {
437: int tag;
438:
439: tag = receive_ox_tag(fd_rw);
440: switch(tag) {
441: case OX_DATA:
442: push(receive_cmo(fd_rw));
443: if (Debug) show_stack_top();
444: break;
445: case OX_COMMAND:
446: if (Debug) show_stack_top();
447: receive_and_execute_sm_command();
448: break;
449: default:
450: ;
451: }
452: return 0;
453: }
454:
1.2 takayama 455: jmp_buf Ox_env;
1.3 takayama 456: int Ox_intr_usr1=0;
1.2 takayama 457: void usr1_handler(int sig)
458: {
1.3 takayama 459: Ox_intr_usr1=1;
460: longjmp(Ox_env,1);
461: }
462: void restart() {
463: Ox_intr_usr1=0;
1.2 takayama 464: longjmp(Ox_env,1);
465: }
466:
1.3 takayama 467: void myhandler(const char *reason,const char *file,int line, int gsl_errno) {
468: cmo *m;
469: FILE *fp;
470: char logname[1024];
471: sprintf(logname,"/tmp/ox_gsl-%d.txt",(int) getpid());
472: fp = fopen(logname,"w");
473: fprintf(fp,"%d\n",gsl_errno);
474: fprintf(fp,"%d\n",line);
475: if (file != NULL) fprintf(fp,"%s\n",file); else fprintf(fp,"file?\n");
476: if (reason != NULL) fprintf(fp,"%s\n",reason); else fprintf(fp,"reason?\n");
1.4 takayama 477: fflush(NULL); fclose(fp);
1.3 takayama 478: // m = make_error2(reason,file,line,gsl_errno);
479: // send_ox_cmo(fd_rw, m); ox_flush(fd_rw);
480: // send error packet even it is not asked. Todo, OK? --> no
481: restart();
482: }
483: void push_error_from_file() {
484: FILE *fp;
485: #define BUF_SIZE 1024
486: char logname[BUF_SIZE];
487: char cmd[BUF_SIZE];
488: char file[BUF_SIZE];
489: char reason[BUF_SIZE];
490: int gsl_errno, line;
491: cmo *m;
492: fprintf(stderr,"push_error_from_file()\n");
493: sprintf(logname,"/tmp/ox_gsl-%d.txt",(int) getpid());
1.4 takayama 494: fp = fopen(logname,"r");
495: if (fp == NULL) {
496: fprintf(stderr,"open %s is failed\n",logname); return;
497: }
1.3 takayama 498: fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&gsl_errno);
499: fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&line);
1.4 takayama 500: #define remove_newline(s) {char *tmp_pos; if ((tmp_pos=strchr(s,'\n')) != NULL) *tmp_pos = '\0';}
501: fgets(file,BUF_SIZE-2,fp); remove_newline(file);
502: fgets(reason,BUF_SIZE-2,fp); remove_newline(reason);
1.3 takayama 503: fclose(fp);
504: m = make_error2(reason,file,line,gsl_errno);
505: push(m);
506: sprintf(cmd,"rm -f %s",logname);
507: system(cmd);
508: }
1.1 takayama 509: int main()
510: {
1.2 takayama 511: if ( setjmp(Ox_env) ) {
1.3 takayama 512: fprintf(stderr,"resetting libgsl ...");
1.2 takayama 513: initialize_stack();
1.3 takayama 514: if (Ox_intr_usr1) {
515: fprintf(stderr,"and sending OX_SYNC_BALL...");
516: send_ox_tag(fd_rw,OX_SYNC_BALL);
517: }
1.2 takayama 518: fprintf(stderr,"done\n");
1.3 takayama 519: Ox_intr_usr1=0;
520: push_error_from_file();
1.2 takayama 521: }else{
1.1 takayama 522: ox_stderr_init(stderr);
523: initialize_stack();
524: init_gc();
525: fd_rw = oxf_open(3);
526: oxf_determine_byteorder_server(fd_rw);
1.2 takayama 527: }
1.10 ohara 528: #if defined(__CYGWIN__)
529: void *mysignal(int sig,void (*handler)(int m));
530: mysignal(SIGUSR1,usr1_handler);
531: #else
1.2 takayama 532: signal(SIGUSR1,usr1_handler);
1.10 ohara 533: #endif
1.2 takayama 534:
535: while(1) {
536: receive();
537: }
538: return(0);
1.1 takayama 539: }
1.13 takayama 540:
541: cmo *element_of_at(cmo *list,int k) {
542: int length;
543: static cmo * saved_list = NULL;
544: static cmo **dic;
545: int i;
546: cell *cellp;
547: if (list == NULL) {
548: ox_printf("element_of_at: list is NULL.\n");
549: return( (cmo *)NULL);
550: }
551: if (list->tag != CMO_LIST) {
552: ox_printf("element_of_at: list is not list.\n");
553: return((cmo *)NULL);
554: }
555: length = list_length((cmo_list *)list);
556: if ((k < 0) || (k >= length)) {
557: ox_printf("element_of_at: out of bound length=%d, k=%d.\n",length,k);
558: return((cmo *)NULL);
559: }
560: if (list == saved_list) return(dic[k]);
561: saved_list = list;
562: dic = (cmo **)GC_malloc(sizeof(cmo *)*(length+1));
563: if (dic == NULL) return((cmo *)NULL); // no more memory.
564: cellp = list_first((cmo_list *)list);
565: for (i=0; i<length; i++) {
566: dic[i] = cellp->cmo;
567: cellp = list_next(cellp);
568: }
569: return(dic[k]);
570: }
1.15 ! takayama 571:
! 572: int get_length(cmo *c) {
! 573: if (c->tag != CMO_LIST) {
! 574: return(-1);
! 575: }
! 576: return(list_length((cmo_list *)c));
! 577: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ox_gsl.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / ox_gsl