[BACK]Return to ox_gsl.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / ox_gsl

Annotation of OpenXM/src/ox_gsl/ox_gsl.c, Revision 1.12

1.12    ! takayama    1: /* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.11 2018/06/06 07:40:32 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: }
                    226:
                    227: void my_add_double() {
                    228:   double x,y;
                    229:   pop();
                    230:   y = get_double();
                    231:   x = get_double();
                    232:   push((cmo *)new_cmo_double(x+y));
                    233: }
                    234:
                    235: double *get_double_list(int *length) {
                    236:   cmo *c;
                    237:   cmo *entry;
                    238:   cell *cellp;
                    239:   double *d;
                    240:   int n,i;
                    241:   c = pop();
                    242:   if (c->tag != CMO_LIST) {
1.3       takayama  243: //    make_error2("get_double_list",NULL,0,-1);
1.2       takayama  244:     *length=-1; return(0);
                    245:   }
                    246:   n = *length = list_length((cmo_list *)c);
                    247:   d = (double *) GC_malloc(sizeof(double)*(*length+1));
                    248:   cellp = list_first((cmo_list *)c);
                    249:   entry = cellp->cmo;
                    250:   for (i=0; i<n; i++) {
                    251:     if (Debug) {
                    252:       printf("entry[%d]=",i); print_cmo(entry); printf("\n");
                    253:     }
                    254:     if (entry->tag == CMO_INT32) {
                    255:       d[i]=( (double) (((cmo_int32 *)entry)->i) );
                    256:     }else if (entry->tag == CMO_IEEE_DOUBLE_FLOAT) {
                    257:       d[i]=((cmo_double *)entry)->d;
                    258:     }else if (entry->tag == CMO_ZZ) {
                    259:       d[i]=( (double) mpz_get_si(((cmo_zz *)entry)->mpz));
                    260:     }else if (entry->tag == CMO_NULL) {
                    261:       d[i]= 0;
                    262:     }else {
                    263:       fprintf(stderr,"entries of the list should be int32 or zz or double\n");
                    264:       *length = -1;
1.3       takayama  265:       myhandler("get_double_list",NULL,0,-1);
1.2       takayama  266:       return(NULL);
                    267:     }
                    268:     cellp = list_next(cellp);
                    269:     entry = cellp->cmo;
                    270:   }
                    271:   return(d);
                    272: }
                    273: void show_double_list() {
                    274:   int n;
                    275:   double *d;
                    276:   int i;
                    277:   pop(); // pop argument number;
                    278:   d = get_double_list(&n);
1.3       takayama  279:   if (n < 0) fprintf(stderr,"Error in the double list\n");
1.2       takayama  280:   printf("show_double_list: length=%d\n",n);
                    281:   for (i=0; i<n; i++) {
                    282:     printf("%lg, ",d[i]);
                    283:   }
                    284:   printf("\n");
                    285: }
                    286:
                    287: char *get_string() {
                    288:   cmo *c;
                    289:   c = pop();
                    290:   if (c->tag == CMO_STRING) {
                    291:     return (((cmo_string *)c)->s);
                    292:   }
1.3       takayama  293:   // make_error2(-1);
1.2       takayama  294:   return(NULL);
                    295: }
1.1       takayama  296:
1.7       takayama  297: void test_ox_eval() {
1.5       takayama  298:   cmo *c;
                    299:   double d=0;
                    300:   pop();
1.7       takayama  301:   c=pop();
1.5       takayama  302:   if (Debug) {
1.7       takayama  303:     ox_printf("cmo *c="); print_cmo(c); ox_printf("\n");
1.5       takayama  304:   }
1.6       ohara     305:   init_dic();
1.5       takayama  306:   register_entry("x",1.25);
1.7       takayama  307:   if (eval_cmo(c,&d) == 0) myhandler("eval_cmo failed",NULL,0,-1);
1.5       takayama  308:   push((cmo *)new_cmo_double(d));
                    309: }
                    310:
1.1       takayama  311: int sm_executeFunction()
                    312: {
                    313:     cmo_string *func = (cmo_string *)pop();
                    314:     if (func->tag != CMO_STRING) {
1.3       takayama  315:         push(make_error2("sm_executeFunction, not CMO_STRING",NULL,0,-1));
1.1       takayama  316:         return -1;
                    317:     }
1.9       takayama  318:     init_dic();
1.2       takayama  319:     // Test functions
1.1       takayama  320:     if (strcmp(func->s, "add_int32") == 0) {
                    321:         my_add_int32();
1.2       takayama  322:     }else if (strcmp(func->s,"add_double")==0) {
                    323:         my_add_double();
                    324:     }else if (strcmp(func->s,"show_double_list")==0) {
                    325:         show_double_list();
1.3       takayama  326:     }else if (strcmp(func->s,"restart")==0) {
                    327:         pop(); restart();
1.5       takayama  328:     }else if (strcmp(func->s,"test_ox_eval")==0) {
                    329:         test_ox_eval();
1.2       takayama  330:     // The following functions are defined in call_gsl.c
                    331:     }else if (strcmp(func->s,"gsl_sf_lngamma_complex_e")==0) {
                    332:         call_gsl_sf_lngamma_complex_e();
1.8       takayama  333:     }else if (strcmp(func->s,"gsl_integration_qags")==0) {
                    334:         call_gsl_integration_qags();
1.11      takayama  335:     }else if (strcmp(func->s,"gsl_monte_plain_integrate")==0) {
1.12    ! takayama  336:         call_gsl_monte_plain_miser_vegas_integrate(0);
        !           337:     }else if (strcmp(func->s,"gsl_monte_miser_integrate")==0) {
        !           338:         call_gsl_monte_plain_miser_vegas_integrate(1);
        !           339:     }else if (strcmp(func->s,"gsl_monte_vegas_integrate")==0) {
        !           340:         call_gsl_monte_plain_miser_vegas_integrate(2);
1.1       takayama  341:     }else {
1.3       takayama  342:         push(make_error2("sm_executeFunction, unknown function",NULL,0,-1));
1.1       takayama  343:         return -1;
1.2       takayama  344:     }
                    345:     return(0);
1.1       takayama  346: }
                    347:
                    348:
                    349: int receive_and_execute_sm_command()
                    350: {
                    351:     int code = receive_int32(fd_rw);
                    352:     switch(code) {
                    353:     case SM_popCMO:
                    354:         sm_popCMO();
                    355:         break;
                    356:     case SM_executeFunction:
                    357:         sm_executeFunction();
                    358:         break;
                    359:     case SM_mathcap:
                    360:         sm_mathcap();
                    361:         break;
                    362:     case SM_setMathCap:
                    363:         pop();
                    364:         break;
                    365:     default:
                    366:                 ;
                    367:     }
1.2       takayama  368:     return(0);
1.1       takayama  369: }
                    370:
                    371: int receive()
                    372: {
                    373:     int tag;
                    374:
                    375:     tag = receive_ox_tag(fd_rw);
                    376:     switch(tag) {
                    377:     case OX_DATA:
                    378:         push(receive_cmo(fd_rw));
                    379:         if (Debug) show_stack_top();
                    380:         break;
                    381:     case OX_COMMAND:
                    382:         if (Debug) show_stack_top();
                    383:         receive_and_execute_sm_command();
                    384:         break;
                    385:     default:
                    386:                 ;
                    387:     }
                    388:     return 0;
                    389: }
                    390:
1.2       takayama  391: jmp_buf Ox_env;
1.3       takayama  392: int Ox_intr_usr1=0;
1.2       takayama  393: void usr1_handler(int sig)
                    394: {
1.3       takayama  395:   Ox_intr_usr1=1;
                    396:   longjmp(Ox_env,1);
                    397: }
                    398: void restart() {
                    399:   Ox_intr_usr1=0;
1.2       takayama  400:   longjmp(Ox_env,1);
                    401: }
                    402:
1.3       takayama  403: void myhandler(const char *reason,const char *file,int line, int gsl_errno) {
                    404:   cmo *m;
                    405:   FILE *fp;
                    406:   char logname[1024];
                    407:   sprintf(logname,"/tmp/ox_gsl-%d.txt",(int) getpid());
                    408:   fp = fopen(logname,"w");
                    409:   fprintf(fp,"%d\n",gsl_errno);
                    410:   fprintf(fp,"%d\n",line);
                    411:   if (file != NULL) fprintf(fp,"%s\n",file); else fprintf(fp,"file?\n");
                    412:   if (reason != NULL) fprintf(fp,"%s\n",reason); else fprintf(fp,"reason?\n");
1.4       takayama  413:   fflush(NULL); fclose(fp);
1.3       takayama  414:   // m = make_error2(reason,file,line,gsl_errno);
                    415:   //  send_ox_cmo(fd_rw, m);  ox_flush(fd_rw);
                    416:   // send error packet even it is not asked. Todo, OK? --> no
                    417:   restart();
                    418: }
                    419: void push_error_from_file() {
                    420:   FILE *fp;
                    421: #define BUF_SIZE 1024
                    422:   char logname[BUF_SIZE];
                    423:   char cmd[BUF_SIZE];
                    424:   char file[BUF_SIZE];
                    425:   char reason[BUF_SIZE];
                    426:   int gsl_errno, line;
                    427:   cmo *m;
                    428:   fprintf(stderr,"push_error_from_file()\n");
                    429:   sprintf(logname,"/tmp/ox_gsl-%d.txt",(int) getpid());
1.4       takayama  430:   fp = fopen(logname,"r");
                    431:   if (fp == NULL) {
                    432:     fprintf(stderr,"open %s is failed\n",logname); return;
                    433:   }
1.3       takayama  434:   fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&gsl_errno);
                    435:   fgets(cmd,BUF_SIZE-2,fp); sscanf(cmd,"%d",&line);
1.4       takayama  436: #define remove_newline(s) {char *tmp_pos; if ((tmp_pos=strchr(s,'\n')) != NULL) *tmp_pos = '\0';}
                    437:   fgets(file,BUF_SIZE-2,fp);  remove_newline(file);
                    438:   fgets(reason,BUF_SIZE-2,fp); remove_newline(reason);
1.3       takayama  439:   fclose(fp);
                    440:   m = make_error2(reason,file,line,gsl_errno);
                    441:   push(m);
                    442:   sprintf(cmd,"rm -f %s",logname);
                    443:   system(cmd);
                    444: }
1.1       takayama  445: int main()
                    446: {
1.2       takayama  447:   if ( setjmp(Ox_env) ) {
1.3       takayama  448:     fprintf(stderr,"resetting libgsl ...");
1.2       takayama  449:     initialize_stack();
1.3       takayama  450:     if (Ox_intr_usr1) {
                    451:       fprintf(stderr,"and sending OX_SYNC_BALL...");
                    452:       send_ox_tag(fd_rw,OX_SYNC_BALL);
                    453:     }
1.2       takayama  454:     fprintf(stderr,"done\n");
1.3       takayama  455:     Ox_intr_usr1=0;
                    456:     push_error_from_file();
1.2       takayama  457:   }else{
1.1       takayama  458:     ox_stderr_init(stderr);
                    459:     initialize_stack();
                    460:     init_gc();
                    461:     fd_rw = oxf_open(3);
                    462:     oxf_determine_byteorder_server(fd_rw);
1.2       takayama  463:   }
1.10      ohara     464: #if defined(__CYGWIN__)
                    465:   void *mysignal(int sig,void (*handler)(int m));
                    466:   mysignal(SIGUSR1,usr1_handler);
                    467: #else
1.2       takayama  468:   signal(SIGUSR1,usr1_handler);
1.10      ohara     469: #endif
1.2       takayama  470:
                    471:   while(1) {
                    472:     receive();
                    473:   }
                    474:   return(0);
1.1       takayama  475: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>