=================================================================== RCS file: /home/cvs/OpenXM/src/ox_gsl/ox_gsl.c,v retrieving revision 1.13 retrieving revision 1.18 diff -u -p -r1.13 -r1.18 --- OpenXM/src/ox_gsl/ox_gsl.c 2018/06/07 11:13:05 1.13 +++ 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.12 2018/06/07 01:53:33 takayama Exp $ +/* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.17 2019/10/23 07:00:43 takayama Exp $ */ #include @@ -10,6 +10,7 @@ #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; @@ -223,6 +224,30 @@ double get_double() 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; @@ -270,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; @@ -338,6 +401,12 @@ int sm_executeFunction() 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; @@ -420,7 +489,7 @@ 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; @@ -503,4 +572,11 @@ cmo *element_of_at(cmo *list,int k) { 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)); }