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