lib/cffi-glue.c

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