Return to ox_gsl.c CVS log

Up to [local] / OpenXM / src / ox_gsl

The complex gamma function implemented in GSL is installed.
Ex.
ox_launch(0,"/Users/nobuki/OX4/OpenXM/src/ox_gsl/ox_gsl");
0
[1916] ox_cmo_rpc(0,"gsl_sf_lngamma_complex_e",3,0); ox_pop_cmo(0);
0
[1917] [0.693147,0,0]   [log( |Gamma(z)| ), arg( Gamma(z)), status]
[1918] ox_cmo_rpc(0,"gsl_sf_lngamma_complex_e",3,2); ox_pop_cmo(0);
0
[1919] [4.44406e+15,2.32392,0]
CVS ----------------------------------------------------------------------

/* $OpenXM: OpenXM/src/ox_gsl/ox_gsl.c,v 1.2 2018/03/29 11:52:18 takayama Exp $
*/

#include <stdio.h>
#include <stdlib.h>
#include <setjmp.h>
#include <string.h>
#include "ox_toolkit.h"

OXFILE *fd_rw;

#define INIT_S_SIZE 2048
#define EXT_S_SIZE  2048

static int stack_size = 0;
static int stack_pointer = 0;
static cmo **stack = NULL;

int Debug=1;

void show_stack_top() {
  cmo *data;
  if (stack_pointer > 0) {
    data=stack[stack_pointer-1];
    print_cmo(data); printf("\n");
  }else {
    printf("The stack is empty.\n");
  }
}

void init_gc()
{
  GC_INIT();
}

void initialize_stack()
{
    stack_pointer = 0;
    stack_size = INIT_S_SIZE;
    stack = malloc(stack_size*sizeof(cmo*));
}

static void extend_stack()
{
    int size2 = stack_size + EXT_S_SIZE;
    cmo **stack2 = malloc(size2*sizeof(cmo*));
    memcpy(stack2, stack, stack_size*sizeof(cmo *));
    free(stack);
    stack = stack2;
    stack_size = size2;
}

void push(cmo* m)
{
    stack[stack_pointer] = m;
    stack_pointer++;
    if (stack_pointer >= stack_size) {
        extend_stack();
    }
}

cmo* pop()
{
    if (stack_pointer > 0) {
        stack_pointer--;
        return stack[stack_pointer];
    }
    return new_cmo_null();
}

void pops(int n)
{
    stack_pointer -= n;
    if (stack_pointer < 0) {
        stack_pointer = 0;
    }
}

#define OX_GSL_VERSION 2018032901
#define ID_STRING  "2018/03/29 13:56:00"

int sm_mathcap()
{
  int available_cmo[]={
    CMO_NULL,
    CMO_INT32,
//    CMO_DATUM,
    CMO_STRING,
    CMO_MATHCAP,
    CMO_LIST,
//    CMO_MONOMIAL32,
    CMO_ZZ,
//    CMO_QQ,
    CMO_BIGFLOAT32,
    CMO_COMPLEX,
    CMO_IEEE_DOUBLE_FLOAT,
    CMO_ZERO,
//    CMO_DMS_GENERIC,
//    CMO_RING_BY_NAME,
//    CMO_INDETERMINATE,
//    CMO_DISTRIBUTED_POLYNOMIAL,
//    CMO_RECURSIVE_POLYNOMIAL,
//    CMO_POLYNOMIAL_IN_ONE_VARIABLE,
    CMO_ERROR2,
    0};
  int available_sm_command[]={
    SM_popCMO,
    SM_popString,
    SM_mathcap,
    SM_pops,
//    SM_executeStringByLocalParser,
    SM_executeFunction,
    SM_setMathCap,
    SM_shutdown,
    SM_control_kill,
    SM_control_reset_connection,
    SM_control_spawn_server,
    SM_control_terminate_server,
    0};
    mathcap_init(OX_GSL_VERSION, ID_STRING, "ox_gsl", available_cmo,available_sm_command);
    push((cmo *)oxf_cmo_mathcap(fd_rw));
    return 0;
}

int sm_popCMO()
{
    cmo* m = pop();

    if (m != NULL) {
        send_ox_cmo(fd_rw, m);
        return 0;
    }
    return SM_popCMO;
}

cmo *make_error2(int code)
{
    fprintf(stderr,"make_error2: not implemented.\n");
    return ((cmo *)new_cmo_int32(-1));
}

int get_i()
{
    cmo *c = pop();
    if (c->tag == CMO_INT32) {
        return ((cmo_int32 *)c)->i;
    }else if (c->tag == CMO_ZZ) {
        return mpz_get_si(((cmo_zz *)c)->mpz);
    }else if (c->tag == CMO_NULL) {
        return(0);
    }else if (c->tag == CMO_ZERO) {
        return(0);
    }
    make_error2(-1);
    return 0;
}

void get_xy(int *x, int *y)
{
    pop();
    *x = get_i();
    *y = get_i();
}

void my_add_int32()
{
    int x, y;
    get_xy(&x, &y);
    push((cmo *)new_cmo_int32(x+y));
}

double get_double()
{
    cmo *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) {
        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);
    }
    make_error2(-1);
    return 0;
}

void my_add_double() {
  double x,y;
  pop();
  y = get_double();
  x = get_double();
  push((cmo *)new_cmo_double(x+y));
}

double *get_double_list(int *length) {
  cmo *c;
  cmo *entry;
  cell *cellp;
  double *d;
  int n,i;
  c = pop();
  if (c->tag != CMO_LIST) {
    make_error2(-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; i<n; i++) {
    if (Debug) {
      printf("entry[%d]=",i); print_cmo(entry); printf("\n");
    }
    if (entry->tag == 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");
      make_error2(-1);
      *length = -1;
      return(NULL);
    }
    cellp = list_next(cellp);
    entry = cellp->cmo;
  }
  return(d);
}
void show_double_list() {
  int n;
  double *d;
  int i;
  pop(); // pop argument number;
  d = get_double_list(&n);
  printf("show_double_list: length=%d\n",n);
  for (i=0; i<n; i++) {
    printf("%lg, ",d[i]);
  }
  printf("\n");
}

char *get_string() {
  cmo *c;
  c = pop();
  if (c->tag == CMO_STRING) {
    return (((cmo_string *)c)->s);
  }
  make_error2(-1);
  return(NULL);
}

int sm_executeFunction()
{
    cmo_string *func = (cmo_string *)pop();
    if (func->tag != CMO_STRING) {
        push(make_error2(0));
        return -1;
    }
    // Test functions
    if (strcmp(func->s, "add_int32") == 0) {
        my_add_int32();
    }else if (strcmp(func->s,"add_double")==0) {
        my_add_double();
    }else if (strcmp(func->s,"show_double_list")==0) {
        show_double_list();
    // 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 {
        push(make_error2(0));
        return -1;
    } 
    return(0);
}


int receive_and_execute_sm_command()
{
    int code = receive_int32(fd_rw);
    switch(code) {
    case SM_popCMO:
        sm_popCMO();
        break;
    case SM_executeFunction:
        sm_executeFunction();
        break;
    case SM_mathcap:
        sm_mathcap();
        break;
    case SM_setMathCap:
        pop();
        break;
    default:
                ;
    }
    return(0);
}

int receive()
{
    int tag;

    tag = receive_ox_tag(fd_rw);
    switch(tag) {
    case OX_DATA:
        push(receive_cmo(fd_rw));
        if (Debug) show_stack_top();
        break;
    case OX_COMMAND:
        if (Debug) show_stack_top();
        receive_and_execute_sm_command();
        break;
    default:
                ;
    }
    return 0;
}

jmp_buf Ox_env;

void usr1_handler(int sig)
{
  longjmp(Ox_env,1);
}

int main()
{
  if ( setjmp(Ox_env) ) {
    fprintf(stderr,"resetting libgsl and sending OX_SYNC_BALL...");
    initialize_stack();
    send_ox_tag(fd_rw,OX_SYNC_BALL);
    fprintf(stderr,"done\n");
  }else{
    ox_stderr_init(stderr);
    initialize_stack();
    init_gc();
    fd_rw = oxf_open(3);
    oxf_determine_byteorder_server(fd_rw);
  }
  signal(SIGUSR1,usr1_handler);
  
  while(1) {
    receive();
  }
  return(0);
}