lib/cffi-glue.c

   1 /*
   2 #ifdef __APPLE__
   3 #include <stdlib.h>
   4 #else
   5 #include <malloc.h>
   6 #endif
   7 */
   8 #include <stdlib.h>
   9 #include <malloc.h>
  10
  11
  12
  13 #include <string.h>
  14
  15 #include "linalg.h"
  16
  17 typedef char *PTR;
  18
  19 extern double unirand(), gamma();
  20 extern double normalcdf(), normalquant(), normaldens(), normalrand();
  21 extern double bnormcdf();
  22 extern double cauchycdf(), cauchyquant(), cauchydens(), cauchyrand();
  23 extern double gammacdf(), gammaquant(), gammadens(), gammarand();
  24 extern double chisqcdf(), chisqquant(), chisqdens(), chisqrand();
  25 extern double betacdf(), betaquant(), betadens(), betarand();
  26 extern double tcdf(), tquant(), tdens(), trand();
  27 extern double fcdf(), fquant(), fdens(), frand();
  28 extern double poissoncdf(), poissonpmf();
  29 extern long poissonquant(), poissonrand();
  30 extern double binomialcdf(), binomialpmf();
  31 extern long binomialquant(), binomialrand();
  32 extern double rcondest_front();
  33
  34 void xlfail(char *);
  35
  36 extern int chol_decomp_front();
  37 extern int lu_decomp_front();
  38 extern int lu_solve_front();
  39 extern int lu_inverse_front();
  40 extern int sv_decomp_front();
  41 extern int qr_decomp_front();
  42 extern int make_rotation_front();
  43 extern int eigen_front();
  44 extern int range_to_rseq();
  45 extern int spline_front();
  46 extern int kernel_dens_front();
  47 extern int kernel_smooth_front();
  48 extern int base_lowess_front();
  49 extern int fft_front();
  50 extern int numgrad_front();
  51 extern int numhess_front();
  52 extern int minfo_maximize();
  53
  54 /***********************************************************************/
  55 /****                        Basic Utilities                        ****/
  56 /***********************************************************************/
  57
  58 /***********************************************************************/
  59 /**                    Callback Support Functions                     **/
  60 /***********************************************************************/
  61
  62 static long ccl_integer_value;
  63 static double ccl_double_value;
  64 static PTR ccl_ptr_value;
  65 void ccl_store_integer(long x)  { ccl_integer_value = x; }
  66 void ccl_store_double(double x)  { ccl_double_value = x; }
  67 void ccl_store_ptr(PTR x) { ccl_ptr_value = x; }
  68
  69 /***************************************************************************/
  70 /**                       Lisp-Managed Calloc/Free                        **/
  71 /***************************************************************************/
  72
  73 #ifdef DODO
  74 static void (*new_ptr)();
  75 static void (*free_ptr)();
  76
  77 register_new_ptr(f) void (*f)(); { new_ptr = f; }
  78 register_free_ptr(f) void (*f)(); { free_ptr = f; }
  79
  80 char* calloc(size_t n, size_t m)
  81 {
  82   size_t i, N = n * m;
  83   char *p;
  84
  85   (*new_ptr)(N);
  86   p = (char *) ccl_ptr_value;
  87   for (i = 0; i < N; i++) p[i] = 0;
  88   return(p);
  89 }
  90
  91 malloc() { xlfail("malloc not available yet"); }
  92 realloc() { xlfail("realloc not available yet"); }
  93
  94 void free(char *p)
  95 {
  96   (*free_ptr)(p);
  97 }
  98 #endif /* DODO*/
  99
 100 /***************************************************************************/
 101 /**                                                                       **/
 102 /**                     Storage Allocation Functions                      **/
 103 /**                                                                       **/
 104 /***************************************************************************/
 105
 106 PTR
 107 la_base_allocate(size_t n, size_t m)
 108 {
 109   char *p = calloc(n, m);
 110   if (p == nil) xlfail("allocation failed");
 111   return((PTR) p);
 112 }
 113
 114 long
 115 la_base_free_alloc(PTR p)
 116 {
 117   if (p) free((char *) p);
 118   return(0);
 119 }
 120
 121 size_t
 122 la_mode_size(int mode)
 123 {
 124   switch (mode) {
 125   case IN: return(sizeof(long));
 126   case RE: return(sizeof(double));
 127   case CX: return(sizeof(Complex));
 128   }
 129   return(0);
 130 }
 131
 132 /***************************************************************************/
 133 /**                                                                       **/
 134 /**                    Callbacks for Internal Storage                     **/
 135 /**                                                                       **/
 136 /***************************************************************************/
 137
 138 int (*ccl_la_allocate)(), (*ccl_la_free_alloc)();
 139
 140 void register_la_allocate(f) int (*f)(); { ccl_la_allocate = f; }
 141 void register_la_free_alloc(f) int (*f)(); { ccl_la_free_alloc = f; }
 142
 143 PTR
 144 la_allocate(size_t n, size_t m)
 145 {
 146   (*ccl_la_allocate)(n, m);
 147   return(ccl_ptr_value);
 148 }
 149
 150 void
 151 la_free_alloc(PTR p)
 152 {
 153   (*ccl_la_free_alloc)(p);
 154 }
 155
 156 /***************************************************************************/
 157 /**                                                                       **/
 158 /**                     Storage Access Functions                          **/
 159 /**                                                                       **/
 160 /***************************************************************************/
 161
 162 long
 163 la_get_integer(PTR p, size_t i)
 164 {
 165   return(*(((long *) p) + i));
 166 }
 167
 168 double
 169 la_get_double(PTR p, size_t i)
 170 {
 171   return(*(((double *) p) + i));
 172 }
 173
 174 double
 175 la_get_complex_real(PTR p, size_t i)
 176 {
 177   Complex *c = ((Complex *) p) + i;
 178   return(c->real);
 179 }
 180
 181 double
 182 la_get_complex_imag(PTR p, size_t i)
 183 {
 184   Complex *c = ((Complex *) p) + i;
 185   return(c->imag);
 186 }
 187
 188 PTR
 189 la_get_pointer(PTR p, size_t i)
 190 {
 191   return(*(((PTR *) p) + i));
 192 }
 193
 194 /***************************************************************************/
 195 /**                                                                       **/
 196 /**                     Storage Mutation Functions                        **/
 197 /**                                                                       **/
 198 /***************************************************************************/
 199
 200 int
 201 la_put_integer(PTR p, size_t i, long x)
 202 {
 203   *(((long *) p) + i) = x;
 204   return(0);
 205 }
 206
 207 int la_put_double(PTR p, size_t i, double x)
 208 {
 209   *(((double *) p) + i) = x;
 210   return(0);
 211 }
 212
 213 int
 214 la_put_complex(PTR p, size_t i, double x, double y)
 215 {
 216   Complex *c = ((Complex *) p) + i;
 217   c->real = x;
 218   c->imag = y;
 219   return(0);
 220 }
 221
 222 int
 223 la_put_pointer(PTR p, size_t i, PTR x)
 224 {
 225   *(((PTR *) p) + i) = x;
 226   return(0);
 227 }
 228
 229 /***********************************************************************/
 230 /**                                                                   **/
 231 /**               XLISP Internal Error Message Emulation              **/
 232 /**                                                                   **/
 233 /***********************************************************************/
 234
 235 char buf[1000];
 236
 237 static int (*ccl_set_buf_char_fptr)();
 238 void register_set_buf_char(f) int (*f)(); { ccl_set_buf_char_fptr = f; }
 239 void set_buf_char(int n, int c)  { (*ccl_set_buf_char_fptr)(n, c); }
 240
 241 static int (*ccl_print_buffer)();
 242 void register_print_buffer(f) int (*f)(); { ccl_print_buffer = f; }
 243 void print_buffer(int n, int m) { (*ccl_print_buffer)(n, m); }
 244
 245 static int bufpos = 0;
 246
 247 static void
 248 resetbuf()
 249 {
 250   bufpos = 0;
 251 }
 252
 253 static void
 254 prbuf(char *s)
 255 {
 256   size_t i, n;
 257
 258   n = strlen(s);
 259   for (i = 0; i <n; i++, bufpos++) set_buf_char(bufpos, s[i]);
 260   set_buf_char(bufpos, 0);
 261 }
 262
 263 void
 264 xlfail(char *s)
 265 {
 266   resetbuf();
 267   prbuf(s);
 268   print_buffer(bufpos, 1);
 269 }
 270
 271 void
 272 stdputstr(char *s)
 273 {
 274   resetbuf();
 275   prbuf(s);
 276   print_buffer(bufpos, 0);
 277 }
 278
 279 void
 280 bufputstr(char *s)
 281 {
 282   resetbuf();
 283   prbuf(s);
 284 }
 285
 286 /***************************************************************************/
 287 /****                                                                   ****/
 288 /*****             Lisp Interfaces to Linear Algebra Routines           ****/
 289 /****                                                                   ****/
 290 /***************************************************************************/
 291
 292 int
 293 ccl_chol_decomp_front(PTR mat, size_t n, PTR dpars)
 294 {
 295   return(chol_decomp_front(mat, n, dpars));
 296 }
 297
 298 int
 299 ccl_lu_decomp_front(PTR mat, size_t n, PTR iv, int mode, PTR dp)
 300 {
 301   return(lu_decomp_front(mat, n, iv, mode, dp));
 302 }
 303
 304 int
 305 ccl_lu_solve_front(PTR a, size_t n, PTR indx, PTR b, int mode)
 306 {
 307   return(lu_solve_front(a, n, indx, b, mode));
 308 }
 309
 310 int
 311 ccl_lu_inverse_front(PTR pmat, size_t n, PTR piv, PTR pv, int mode, PTR pinv)
 312 {
 313   return(lu_inverse_front(pmat, n, piv, pv, mode, pinv));
 314 }
 315
 316 int
 317 ccl_sv_decomp_front(PTR mat, size_t m, size_t n, PTR w, PTR v)
 318 {
 319   return(sv_decomp_front(mat, m, n, w, v));
 320 }
 321
 322
 323 int
 324 ccl_qr_decomp_front(PTR mat, size_t m, size_t n, PTR v, PTR jpvt, int pivot)
 325 {
 326   return(qr_decomp_front(mat, m, n, v, jpvt, pivot));
 327 }
 328
 329 double
 330 ccl_rcondest_front(PTR mat, size_t n)
 331 {
 332   return(rcondest_front(mat, n));
 333 }
 334
 335 int
 336 ccl_make_rotation_front(size_t n, PTR rot, PTR x, PTR y, int use_alpha, double alpha)
 337 {
 338   return(make_rotation_front(n, rot, x, y, use_alpha, alpha));
 339 }
 340
 341 int
 342 ccl_eigen_front(PTR a, size_t n, PTR w, PTR z, PTR fv1)
 343 {
 344   return(eigen_front(a, n, w, z, fv1));
 345 }
 346
 347 int
 348 ccl_range_to_rseq(size_t n, PTR px, size_t ns, PTR pxs)
 349 {
 350   return(range_to_rseq(n, px, ns, pxs));
 351 }
 352
 353 int
 354 ccl_spline_front(size_t n, PTR x, PTR y, size_t ns, PTR xs, PTR ys, PTR work)
 355 {
 356   return(spline_front(n, x, y, ns, xs, ys, work));
 357 }
 358
 359 int
 360 ccl_kernel_dens_front(PTR x, size_t n, double width, PTR xs, PTR ys, size_t ns, int code)
 361 {
 362   return(kernel_dens_front(x, n, width, xs, ys, ns, code));
 363 }
 364
 365 int
 366 ccl_kernel_smooth_front(PTR x, PTR y, size_t n, double width, PTR xs, PTR ys, size_t ns, int code)
 367 {
 368   return(kernel_smooth_front(x, y, n, width, xs, ys, ns, code));
 369 }
 370
 371 int
 372 ccl_base_lowess_front(PTR x, PTR y, size_t n, double f,
 373                       size_t nsteps, double delta, PTR ys, PTR rw, PTR res)
 374 {
 375   return(base_lowess_front(x, y, n, f, nsteps, delta, ys, rw, res));
 376 }
 377
 378 int
 379 ccl_fft_front(size_t n, PTR x, PTR work, int isign)
 380 {
 381   return(fft_front(n, x, work, isign));
 382 }
 383
 384 static int (*ccl_maximize_callback)();
 385
 386 void
 387 register_maximize_callback(f)
 388      int (*f)();
 389 {
 390   ccl_maximize_callback = f;
 391 }
 392
 393 void
 394 maximize_callback(size_t n, PTR px, PTR pfval, PTR pgrad, PTR phess, PTR pderivs)
 395 {
 396   (*ccl_maximize_callback)(n, px, pfval, pgrad, phess, pderivs);
 397 }
 398
 399 void
 400 ccl_numgrad_front(size_t n, PTR px, PTR pgrad, double h, PTR pscale)
 401 {
 402   numgrad_front(n, px, pgrad, h, pscale);
 403 }
 404
 405 void
 406 ccl_numhess_front(size_t n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale)
 407 {
 408   numhess_front(n, px, pf, pgrad, phess, h, pscale);
 409 }
 410
 411 void
 412 ccl_minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose)
 413 {
 414   minfo_maximize(px, pfvals, pscale, pip, pdp, verbose);
 415 }
 416
 417 /***********************************************************************/
 418 /***********************************************************************/
 419 /****                                                               ****/
 420 /****                   Probability Distributions                   ****/
 421 /****                                                               ****/
 422 /***********************************************************************/
 423 /***********************************************************************/
 424
 425 static int (*ccl_uni_fptr)();
 426 void register_uni(f) int (*f)(); { ccl_uni_fptr = f; }
 427 double uni() { (*ccl_uni_fptr)(); return(ccl_double_value); }
 428
 429 double ccl_gamma(double x) { return(gamma(x)); }
 430
 431 double ccl_normalcdf(double x)  { return(normalcdf(x)); }
 432 double ccl_normalquant(double x) { return(normalquant(x)); }
 433 double ccl_normaldens(double x) { return(normaldens(x)); }
 434 double ccl_normalrand() { return(normalrand()); }
 435 double ccl_bnormcdf(double x, double y, double z)  {
 436   return(bnormcdf(x,y,z));
 437 }
 438
 439 double ccl_cauchycdf(double x)  { return(cauchycdf(x)); }
 440 double ccl_cauchyquant(double x)  { return(cauchyquant(x)); }
 441 double ccl_cauchydens(double x)  { return(cauchydens(x)); }
 442 double ccl_cauchyrand() { return(cauchyrand()); }
 443
 444 double ccl_gammacdf(double x, double y)  { return(gammacdf(x, y)); }
 445 double ccl_gammaquant(double x, double y)  { return(gammaquant(x, y)); }
 446 double ccl_gammadens(double x, double y)  { return(gammadens(x, y)); }
 447 double ccl_gammarand(double x)  { return(gammarand(x)); }
 448
 449 double ccl_chisqcdf(double x, double y)  { return(chisqcdf(x, y)); }
 450 double ccl_chisqquant(double x, double y)  { return(chisqquant(x, y)); }
 451 double ccl_chisqdens(double x, double y)  { return(chisqdens(x, y)); }
 452 double ccl_chisqrand(double x)  { return(chisqrand(x)); }
 453
 454 double ccl_betacdf(double x, double y, double z)  { return(betacdf(x, y, z)); }
 455 double ccl_betaquant(double x, double y, double z)  { return(betaquant(x, y, z)); }
 456 double ccl_betadens(double x, double y, double z)  { return(betadens(x, y, z)); }
 457 double ccl_betarand(double x, double y)  { return(betarand(x, y)); }
 458
 459 double ccl_tcdf(double x, double y)  { return(tcdf(x, y)); }
 460 double ccl_tquant(double x, double y)  { return(tquant(x, y)); }
 461 double ccl_tdens(double x, double y)  { return(tdens(x, y)); }
 462 double ccl_trand(double x)  { return(trand(x)); }
 463
 464 double ccl_fcdf(double x, double y, double z)  { return(fcdf(x, y, z)); }
 465 double ccl_fquant(double x, double y, double z)  { return(fquant(x, y, z)); }
 466 double ccl_fdens(double x, double y, double z)  { return(fdens(x, y, z)); }
 467 double ccl_frand(double x, double y)  { return(frand(x, y)); }
 468
 469 double ccl_poissoncdf(double x, double y)  { return(poissoncdf(x, y)); }
 470 long ccl_poissonquant(double x, double y)  { return(poissonquant(x, y)); }
 471 double ccl_poissonpmf(long x, double y)  { return(poissonpmf(x, y)); }
 472 long ccl_poissonrand(double x)  { return(poissonrand(x)); }
 473
 474 double ccl_binomialcdf(double x, long y, double z) { return(binomialcdf(x, y, z)); }
 475 long ccl_binomialquant(double x, long y, double z) { return(binomialquant(x, y, z)); }
 476 double ccl_binomialpmf(long x, long y, double z)  { return(binomialpmf(x, y, z)); }
 477 long ccl_binomialrand(long x, double y)  { return(binomialrand(x, y)); }