=================================================================== RCS file: /home/cvs/OpenXM/src/ox_gsl/ox_gsl.c,v retrieving revision 1.3 retrieving revision 1.18 diff -u -p -r1.3 -r1.18 --- OpenXM/src/ox_gsl/ox_gsl.c 2018/03/30 04:43:16 1.3 +++ OpenXM/src/ox_gsl/ox_gsl.c 2021/03/10 06:51:57 1.18 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.2 2018/03/29 11:52:18 takayama Exp $ +/* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.17 2019/10/23 07:00:43 takayama Exp $ */ #include @@ -6,9 +6,11 @@ #include #include #include +#include #include #include "ox_gsl.h" #include "call_gsl.h" // need only when you bind call_gsl functions. +#include "call_gsl_sf.h" OXFILE *fd_rw; @@ -40,6 +42,7 @@ void gc_free(void *p,size_t size) void init_gc() { GC_INIT(); mp_set_memory_functions(GC_malloc,gc_realloc,gc_free); + init_dic(); // initialize ox_eval.c } void initialize_stack() @@ -110,6 +113,7 @@ int sm_mathcap() // CMO_DISTRIBUTED_POLYNOMIAL, // CMO_RECURSIVE_POLYNOMIAL, // CMO_POLYNOMIAL_IN_ONE_VARIABLE, + CMO_TREE, CMO_ERROR2, 0}; int available_sm_command[]={ @@ -154,11 +158,11 @@ cmo *make_error2(const char *reason,const char *fname, n = 5; argv = (cmo **) GC_malloc(sizeof(cmo *)*n); ms = (cmo *)new_cmo_string("Error"); argv[0] = ms; - if (reason != NULL) s = (char *)GC_malloc(strlen(reason)+1); - else strcpy(s,""); + if (reason != NULL) {s = (char *)GC_malloc(strlen(reason)+1); strcpy(s,reason); + }else strcpy(s,""); ms = (cmo *) new_cmo_string(s); argv[1] = ms; - if (reason != NULL) s = (char *)GC_malloc(strlen(fname)+1); - else strcpy(s,""); + if (fname != NULL) {s = (char *)GC_malloc(strlen(fname)+1); strcpy(s,fname); + }else strcpy(s,""); ms = (cmo *) new_cmo_string(s); argv[2] = ms; err = (cmo *)new_cmo_int32(line); argv[3] = err; err = (cmo *)new_cmo_int32(code); argv[4] = err; @@ -179,7 +183,7 @@ int get_i() }else if (c->tag == CMO_ZERO) { return(0); } - myhandler("get_i, not an integer",NULL,0,-1); + myhandler("get_i: not an integer",NULL,0,-1); return 0; } @@ -208,7 +212,7 @@ double get_double() }else if (c->tag == CMO_ZZ) { if ((mpz_cmp_si(mympz(c),(long int) 0x7fffffff)>0) || (mpz_cmp_si(mympz(c),(long int) -0x7fffffff)<0)) { - myhandler("get_double, out of int32",NULL,0,-1); + myhandler("get_double: out of int32",NULL,0,-1); return(NAN); } return( (double) mpz_get_si(((cmo_zz *)c)->mpz)); @@ -217,9 +221,33 @@ double get_double() }else if (c->tag == CMO_ZERO) { return(0); } - myhandler("get_double, not a double",NULL,0,-1); + myhandler("get_double: not a double",NULL,0,-1); return(NAN); } +/* get_double() will be obsolted and will be replaced by cmo2double(c) */ +double cmo2double(cmo *c) +{ +#define mympz(c) (((cmo_zz *)c)->mpz) + if (c == NULL) c = pop(); + if (c->tag == CMO_INT32) { + return( (double) (((cmo_int32 *)c)->i) ); + }else if (c->tag == CMO_IEEE_DOUBLE_FLOAT) { + return (((cmo_double *)c)->d); // see ox_toolkit.h + }else if (c->tag == CMO_ZZ) { + if ((mpz_cmp_si(mympz(c),(long int) 0x7fffffff)>0) || + (mpz_cmp_si(mympz(c),(long int) -0x7fffffff)<0)) { + myhandler("get_double: out of int32",NULL,0,-1); + return(NAN); + } + return( (double) mpz_get_si(((cmo_zz *)c)->mpz)); + }else if (c->tag == CMO_NULL) { + return(0); + }else if (c->tag == CMO_ZERO) { + return(0); + } + myhandler("cmo2double: not a double",NULL,0,-1); + return(NAN); +} void my_add_double() { double x,y; @@ -267,6 +295,44 @@ double *get_double_list(int *length) { } return(d); } +/* get_double_list will be obsolted and will be replaced by cmo2double_list() */ +double *cmo2double_list(int *length,cmo *c) { + cmo *entry; + cell *cellp; + double *d; + int n,i; + if (c == NULL) c = pop(); + if (c->tag != CMO_LIST) { +// make_error2("get_double_list",NULL,0,-1); + *length=-1; return(0); + } + n = *length = list_length((cmo_list *)c); + d = (double *) GC_malloc(sizeof(double)*(*length+1)); + cellp = list_first((cmo_list *)c); + entry = cellp->cmo; + for (i=0; itag == CMO_INT32) { + d[i]=( (double) (((cmo_int32 *)entry)->i) ); + }else if (entry->tag == CMO_IEEE_DOUBLE_FLOAT) { + d[i]=((cmo_double *)entry)->d; + }else if (entry->tag == CMO_ZZ) { + d[i]=( (double) mpz_get_si(((cmo_zz *)entry)->mpz)); + }else if (entry->tag == CMO_NULL) { + d[i]= 0; + }else { + fprintf(stderr,"entries of the list should be int32 or zz or double\n"); + *length = -1; + myhandler("get_double_list",NULL,0,-1); + return(NULL); + } + cellp = list_next(cellp); + entry = cellp->cmo; + } + return(d); +} void show_double_list() { int n; double *d; @@ -291,6 +357,20 @@ char *get_string() { return(NULL); } +void test_ox_eval() { + cmo *c; + double d=0; + pop(); + c=pop(); + if (Debug) { + ox_printf("cmo *c="); print_cmo(c); ox_printf("\n"); + } + init_dic(); + register_entry("x",1.25); + if (eval_cmo(c,&d) == 0) myhandler("eval_cmo failed",NULL,0,-1); + push((cmo *)new_cmo_double(d)); +} + int sm_executeFunction() { cmo_string *func = (cmo_string *)pop(); @@ -298,6 +378,7 @@ int sm_executeFunction() push(make_error2("sm_executeFunction, not CMO_STRING",NULL,0,-1)); return -1; } + init_dic(); // Test functions if (strcmp(func->s, "add_int32") == 0) { my_add_int32(); @@ -307,9 +388,25 @@ int sm_executeFunction() show_double_list(); }else if (strcmp(func->s,"restart")==0) { pop(); restart(); + }else if (strcmp(func->s,"test_ox_eval")==0) { + test_ox_eval(); // The following functions are defined in call_gsl.c }else if (strcmp(func->s,"gsl_sf_lngamma_complex_e")==0) { call_gsl_sf_lngamma_complex_e(); + }else if (strcmp(func->s,"gsl_integration_qags")==0) { + call_gsl_integration_qags(); + }else if (strcmp(func->s,"gsl_monte_plain_integrate")==0) { + call_gsl_monte_plain_miser_vegas_integrate(0); + }else if (strcmp(func->s,"gsl_monte_miser_integrate")==0) { + call_gsl_monte_plain_miser_vegas_integrate(1); + }else if (strcmp(func->s,"gsl_monte_vegas_integrate")==0) { + call_gsl_monte_plain_miser_vegas_integrate(2); + }else if (strcmp(func->s,"gsl_odeiv_step_rk4")==0) { + call_gsl_odeiv_step("rk4"); + }else if (strcmp(func->s,"gsl_sf_gamma_inc")==0) { + call_gsl_sf_gamma_inc(); + }else if (strcmp(func->s,"gsl_eigen_nonsymmv")==0) { + call_gsl_eigen_nonsymmv(); }else { push(make_error2("sm_executeFunction, unknown function",NULL,0,-1)); return -1; @@ -382,30 +479,32 @@ void myhandler(const char *reason,const char *file,int fprintf(fp,"%d\n",line); if (file != NULL) fprintf(fp,"%s\n",file); else fprintf(fp,"file?\n"); if (reason != NULL) fprintf(fp,"%s\n",reason); else fprintf(fp,"reason?\n"); - fflush(NULL); fclose(fp); // BUG. the contents are deleted after it is closed. + fflush(NULL); fclose(fp); // m = make_error2(reason,file,line,gsl_errno); // send_ox_cmo(fd_rw, m); ox_flush(fd_rw); // send error packet even it is not asked. Todo, OK? --> no restart(); } - void push_error_from_file() { FILE *fp; #define BUF_SIZE 1024 char logname[BUF_SIZE]; - char cmd[BUF_SIZE]; + char cmd[BUF_SIZE+100]; char file[BUF_SIZE]; char reason[BUF_SIZE]; int gsl_errno, line; cmo *m; fprintf(stderr,"push_error_from_file()\n"); sprintf(logname,"/tmp/ox_gsl-%d.txt",(int) getpid()); - fp = fopen(logname,"w"); - if (fp == NULL) return; + fp = fopen(logname,"r"); + if (fp == NULL) { + fprintf(stderr,"open %s is failed\n",logname); return; + } fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&gsl_errno); fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&line); - fgets(file,BUF_SIZE-2,fp); - fgets(reason,BUF_SIZE-2,fp); +#define remove_newline(s) {char *tmp_pos; if ((tmp_pos=strchr(s,'\n')) != NULL) *tmp_pos = '\0';} + fgets(file,BUF_SIZE-2,fp); remove_newline(file); + fgets(reason,BUF_SIZE-2,fp); remove_newline(reason); fclose(fp); m = make_error2(reason,file,line,gsl_errno); push(m); @@ -431,10 +530,53 @@ int main() fd_rw = oxf_open(3); oxf_determine_byteorder_server(fd_rw); } +#if defined(__CYGWIN__) + void *mysignal(int sig,void (*handler)(int m)); + mysignal(SIGUSR1,usr1_handler); +#else signal(SIGUSR1,usr1_handler); +#endif while(1) { receive(); } return(0); +} + +cmo *element_of_at(cmo *list,int k) { + int length; + static cmo * saved_list = NULL; + static cmo **dic; + int i; + cell *cellp; + if (list == NULL) { + ox_printf("element_of_at: list is NULL.\n"); + return( (cmo *)NULL); + } + if (list->tag != CMO_LIST) { + ox_printf("element_of_at: list is not list.\n"); + return((cmo *)NULL); + } + length = list_length((cmo_list *)list); + if ((k < 0) || (k >= length)) { + ox_printf("element_of_at: out of bound length=%d, k=%d.\n",length,k); + return((cmo *)NULL); + } + if (list == saved_list) return(dic[k]); + saved_list = list; + dic = (cmo **)GC_malloc(sizeof(cmo *)*(length+1)); + if (dic == NULL) return((cmo *)NULL); // no more memory. + cellp = list_first((cmo_list *)list); + for (i=0; icmo; + cellp = list_next(cellp); + } + return(dic[k]); +} + +int get_length(cmo *c) { + if (c->tag != CMO_LIST) { + return(-1); + } + return(list_length((cmo_list *)c)); }