wrapper for gsl_ran

[Klaus Schilling <Klaus.Schilling@home.ivm.de>]

I added now support for the ramdom state accessors (more advanced stuff).
For that I forged a new smobtype.
/*--------------------------------------------------------------------*/
/*  GPL and stuff */

#include <stdlib.h>
#include <libguile.h>
#include <guile/gh.h>
#include <gsl_ran.h>
#include <gsl_randist.h>
#include <gsl_ran_switch.h>
#include <libguile/dynl.h>


static ulong ranstate_tag;
static scm_sizet ranstate_die (SCM);
static void *scm2ranstate (SCM);
static void ranstate_type_init (void);

static SCM s_use_taus (void);
static SCM s_use_mrg (void);
static SCM s_use_cmrg (void);
static SCM s_use_uni (void);
static SCM s_use_uni32 (void);
static SCM s_use_zuf (void);
static SCM s_use_rand (void);
static SCM s_seed (SCM);
static SCM s_random (void);
static SCM s_max (void);
static SCM s_uniform ();
static SCM s_gaussian ();
static SCM s_poisson (SCM);
static SCM s_exponential (SCM);
static SCM s_gamma (SCM);
static SCM s_poisson_array (SCM, SCM);
static SCM s_getRandomState (void);
static SCM s_setRandomState (SCM);
static SCM s_seed_wstate (SCM, SCM);
static SCM s_random_wstate (SCM);


scm_sizet
ranstate_die (SCM obj)
{
  void * ptr;
  ptr = (void *) SCM_CDR (obj);
  gh_defer_ints ();
  cfree (ptr);
  gh_allow_ints ();
  return 0;
}

static scm_smobfuns ranstate_funs = {0, ranstate_die, 0, 0};

void *
scm2ranstate (SCM obj)
{ 
  SCM_ASSERT (SCM_NIMP (obj) && SCM_CAR (obj) == ranstate_tag, obj, SCM_ARG3, "scm2ranstate");
  return (void *) SCM_CDR (obj);
}

void
ranstate_type_init ()
{
  ranstate_tag = scm_newsmob (&ranstate_funs);
}



SCM
s_use_taus ()
{
  gh_defer_ints ();
  gsl_ran_use_taus ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_mrg ()
{
  gh_defer_ints ();
  gsl_ran_use_mrg ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_cmrg ()
{
  gh_defer_ints ();
  gsl_ran_use_cmrg ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_uni ()
{
  gh_defer_ints ();
  gsl_ran_use_uni ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_uni32 ()
{
  gh_defer_ints ();
  gsl_ran_use_uni32 ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_rand ()
{
  gh_defer_ints ();
  gsl_ran_use_rand ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_use_zuf ()
{
  gh_defer_ints ();
  gsl_ran_use_zuf ();
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_seed (SCM s_n)
{
  int n;
  n = gh_scm2int (s_n);
  gh_defer_ints ();
  gsl_ran_seed (n);
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM
s_random ()
{
  ulong result;
  gh_defer_ints ();
  result = gsl_ran_random ();
  gh_allow_ints ();
  return gh_ulong2scm (result);
}

SCM
s_max ()
{
  ulong result;
  gh_defer_ints ();
  result = gsl_ran_max ();
  gh_allow_ints ();
  return gh_ulong2scm (result);
}

SCM
s_uniform ()
{
  double result;
  gh_defer_ints ();
  result = gsl_ran_uniform ();
  gh_allow_ints ();
  return gh_double2scm (result);
}

SCM
s_gaussian ()
{
  double result;
  gh_defer_ints ();
  result = gsl_ran_gaussian ();
  gh_allow_ints ();
  return gh_double2scm (result);
}

SCM
s_poisson (SCM s_m)
{
  double m;
  int result;
  m = gh_scm2double (s_m);
  gh_defer_ints ();
  result = gsl_ran_poisson (m);
  gh_allow_ints ();
  return gh_int2scm (result);
}

SCM
s_exponential (SCM s_m)
{
  double m;
  double result;
  m = gh_scm2double (s_m);
  gh_defer_ints ();
  result = gsl_ran_exponential (m);
  gh_allow_ints ();
  return gh_double2scm (result);
}

SCM
s_gamma (SCM s_m)
{
  double m;
  int result;
  m = gh_scm2double (s_m);
  gh_defer_ints ();
  result = gsl_ran_gamma (m);
  gh_allow_ints ();
  return gh_int2scm (result);
}

SCM
s_poisson_array (SCM s_m, SCM s_n)
{
  int n; 
  int *i;
  SCM result;
  double m;
  n = gh_scm2int (s_n);
  m = gh_scm2double (s_m);
  i = (int *) calloc (n, sizeof (int));
  if (i == (int *)NULL) scm_throw (gh_symbol2scm ("memory-exhausted"), SCM_EOL); 
  gsl_ran_poisson_array (m, n, i);
  result = gh_ints2scm (i, n);
  cfree (i);
  return result;   
}

SCM
s_getRandomState (void)
{
  void * ptr;
  SCM scm;
  SCM_NEWCELL (scm);
  gh_defer_ints ();
  ptr = gsl_ran_getRandomState ();
  gh_allow_ints ();
  SCM_SETCAR (scm, ranstate_tag);
  SCM_SETCDR (scm, ptr);
  return scm;
}

SCM 
s_setRandomState (SCM obj)
{
  void *ptr;
  ptr = scm2ranstate (obj);
  gh_defer_ints ();
  gsl_ran_setRandomState (ptr);
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM 
s_seed_wstate (SCM obj, SCM s_i)
{
  int i;
  void *ptr;
  ptr = scm2ranstate (obj);
  i = gh_scm2int (s_i);
  gh_defer_ints ();
  gsl_ran_seed_wstate (ptr, i);
  gh_allow_ints ();
  return SCM_UNSPECIFIED;
}

SCM 
s_random_wstate (SCM obj)
{
  ulong n;
  void *ptr;
  ptr = scm2ranstate (obj);
  gh_defer_ints ();
  n = gsl_ran_random_wstate (ptr);
  gh_allow_ints ();
  return gh_ulong2scm (n);
}






void init_gsl_ran ()
{
  ranstate_type_init ();
  gh_new_procedure ("use-taus", s_use_taus, 0, 0, 0);
  gh_new_procedure ("use-mrg", s_use_mrg, 0, 0, 0);
  gh_new_procedure ("use-cmrg", s_use_cmrg, 0, 0, 0);
  gh_new_procedure ("use-uni", s_use_uni, 0, 0, 0);
  gh_new_procedure ("use-uni32", s_use_uni32, 0, 0, 0);
  gh_new_procedure ("use-rand", s_use_rand, 0, 0, 0);
  gh_new_procedure ("use-zuf", s_use_zuf, 0, 0, 0);
  gh_new_procedure ("seed", s_seed, 1, 0, 0);
  gh_new_procedure ("random", s_random, 0, 0, 0);
  gh_new_procedure ("max", s_max, 0, 0, 0);
  gh_new_procedure ("uniform", s_uniform, 0, 0, 0);
  gh_new_procedure ("gaussian", s_gaussian, 0, 0, 0);
  gh_new_procedure ("poisson", s_poisson, 1, 0, 0);
  gh_new_procedure ("exponential", s_exponential, 1, 0, 0);
  gh_new_procedure ("gamma", s_gamma, 1, 0, 0);
  gh_new_procedure ("poisson-array", s_poisson_array, 2, 0, 0);
  gh_new_procedure ("get-random-state", s_getRandomState, 0, 0, 0);
  gh_new_procedure ("set-random-state", s_setRandomState, 1, 0, 0);
  gh_new_procedure ("seed-wstate", s_seed_wstate, 2, 0, 0);
  gh_new_procedure ("random-wstate", s_random_wstate, 1, 0, 0);
}



void
scm_init_gsl_ran_module ()
{
  scm_register_module_xxx ("gsl ran", init_gsl_ran);
}

/*---------------------------------------------------------------------*/

It is meant to provide a shared extension module gls/ran , from somewhere
in the guile load path, when linked with -lgslrandist and -lgslrandom.


Klaus Schilling