/* vct support 
 *
 * a vct is a Guile "smob" containing a float array and its size
 * we use these in Snd because Guile's floating-point operations are incredibly slow
 * and Snd/CLM are applications where no user is willing to wait on a multiply.
 *
 * C side:
 *   void init_vct(void)                   called to declare the various functions and the vct type in Guile
 *   int vct_p(SCM obj)                    is obj a vct
 *   SCM make_vct(int len, float *data)    make a new vct
 *   vct *get_vct(SCM arg)                 given SCM arg, return vct object
 *   void set_vct_print_length(int val)    set vct print length (default 10)
 *
 * Scheme side:
 *   (make-vct len)                        make new vct
 *   (vct? obj)                            is obj a vct
 *   (vct-ref v index)                     return vct->data[index]
 *   (vct-set! v index val)                vct->data[index] = val
 *   (vct-copy v)                          return a copy of v
 *   (vct-length v)                        return length of vct->data
 *   (vct-add! v1 v2)                      add (element-wise) v2 to v1 (changing v1), return v1
 *   (vct-offset! v1 scl)                  add scl to each element of v1 (changing v1), return v1
 *   (vct-multiply! v1 v2)                 multiply (element-wise, a windowing operation) v1 by v2, changing v1, return v1
 *   (vct-scale! v1 scl)                   multiply each element of v1 by scl, chaning v1, return v1
 *   (vct-fill! v1 val)                    set each element of v1 to val, returning v1
 *   (list->vct lst)                       return vct with elements of lst
 *
 * The intended use is a sort of latter-day array-processing system that handles huge
 * one-dimensional vectors -- fft's, etc.  Some of these functions can be found in
 * snd-gh.c in the Snd package; others can be found in the CLM package, (clm2scm.c).
 */

#if defined(HAVE_CONFIG_H)
  #include "config.h"
#endif
#ifdef DEBUG_MEMORY
  #include <stdlib.h>
  #include "sndlib.h"
#endif

#if HAVE_GUILE

#include "vct.h"

#ifndef CALLOC
  #define CALLOC(a,b)  calloc(a,b)
  #define MALLOC(a)    malloc(a)
  #define FREE(a)      free(a)
  #define REALLOC(a,b) realloc(a,b)
#endif

#define VCT_PRINT_LENGTH 10
#define ERRN1(a,b) SCM_ASSERT((gh_number_p(a)),a,SCM_ARG1,b)
#define ERRN2(a,b) SCM_ASSERT((gh_number_p(a)),a,SCM_ARG2,b)
#define ERRN3(a,b) SCM_ASSERT((gh_number_p(a)),a,SCM_ARG3,b)

#define RTNBOOL(a) return((a) ? SCM_BOOL_T : SCM_BOOL_F)
#define RTNINT(a) return(gh_int2scm(a))
#define RTNFLT(a) return(gh_double2scm(a))

#define GH_TYPE_OF(a) (SCM_TYP16(a))

#ifndef MIN
  #define MIN(a,b) ((a > b) ? (b) : (a))
#endif

static int vct_tag = 0;
static int vct_print_length = VCT_PRINT_LENGTH;
void set_vct_print_length(int len) {vct_print_length = len;}

static SCM mark_vct(SCM obj)
{
  SCM_SETGC8MARK(obj);
  return(SCM_BOOL_F);
}

int vct_p(SCM obj)
{
  return((SCM_NIMP(obj)) && (GH_TYPE_OF(obj) == (SCM)vct_tag));
}

static SCM g_vct_p(SCM obj) {RTNBOOL(vct_p(obj));}

vct *get_vct(SCM arg)
{
  if (vct_p(arg))
    return((vct *)gh_cdr(arg));
  return(NULL);
}

static scm_sizet free_vct(SCM obj)
{
  vct *v = (vct *)gh_cdr(obj);
  if (v->data) FREE(v->data);
  v->data = NULL;
  FREE(v);
  return(0);
}

static int print_vct(SCM obj, SCM port, scm_print_state *pstate)
{
  int len,i;
  char *buf;
  vct *v = (vct *)gh_cdr(obj);
  scm_puts("#<vct",port);
  len = vct_print_length;
  if (len > v->length) len = v->length;
  if (len > 0)
    {
      buf = (char *)CALLOC(32,sizeof(char));
      for (i=0;i<len;i++)
	{
	  sprintf(buf," %.3f",v->data[i]);
	  scm_puts(buf,port);
	}
      if (v->length > vct_print_length)
	scm_puts(" ...",port);
      FREE(buf);
    }
  scm_puts(">",port);
  return(1);
}

static SCM equalp_vct(SCM obj1, SCM obj2)
{
  vct *v1,*v2;
  int i;
  v1 = (vct *)gh_cdr(obj1);
  v2 = (vct *)gh_cdr(obj2);
  if (v1->length != v2->length) return(SCM_BOOL_F);
  for (i=0;i<v1->length;i++)
    if (v1->data[i] != v2->data[i])
      return(SCM_BOOL_F);
  return(SCM_BOOL_T);
}

SCM make_vct(int len, float *data)
{
  SCM ans;
  vct *new_vct;
  new_vct = (vct *)CALLOC(1,sizeof(vct));
  new_vct->length = len;
  new_vct->data = data;
  SCM_NEWCELL(ans);
  SCM_SETCAR(ans,vct_tag);
  SCM_SETCDR(ans,(SCM)new_vct);
  return(ans);
}

#if (!HAVE_MAKE_SMOB_TYPE)
static scm_smobfuns vct_smobfuns = {
  &mark_vct,
  &free_vct,
  &print_vct,
  &equalp_vct};
#endif

static SCM g_make_vct(SCM len)
{
  int size;
  ERRN1(len,S_make_vct);
  size = gh_scm2int(len);
  return(make_vct(size,(float *)CALLOC(size,sizeof(float))));
}

static SCM copy_vct(SCM obj)
{
  vct *v;
  float *copied_data;
  int len,i;
  ERRVCT1(obj,S_vct_copy);
  v = get_vct(obj);
  if (v)
    {
      len = v->length;
      copied_data = (float *)CALLOC(len,sizeof(float));
      for (i=0;i<len;i++) copied_data[i] = v->data[i];
      return(make_vct(len,copied_data));
    }
  return(SCM_BOOL_F);
}

static SCM vct_length(SCM obj)
{
  vct *v = get_vct(obj);
  ERRVCT1(obj,S_vct_length);
  if (v)
    RTNINT(v->length);
  RTNINT(0);
}

static SCM vct_ref(SCM obj, SCM pos)
{
  vct *v = get_vct(obj);
  int loc;
  ERRVCT1(obj,S_vct_ref);
  ERRN2(pos,S_vct_ref);
  if (v)
    {
      loc = gh_scm2int(pos);
      if ((loc >= 0) && (loc < v->length))
	RTNFLT(v->data[loc]);
      else scm_misc_error(S_vct_ref,"invalid index",SCM_LIST2(obj,pos));
    }
  else scm_misc_error(S_vct_ref,"nil vct?",SCM_EOL);
  RTNFLT(0.0);
}

static SCM vct_set(SCM obj, SCM pos, SCM val)
{
  vct *v = get_vct(obj);
  int loc;
  ERRVCT1(obj,S_vct_setB);
  ERRN2(pos,S_vct_setB);
  ERRN3(val,S_vct_setB);
  if (v)
    {
      loc = gh_scm2int(pos);
      if ((loc >= 0) && (loc < v->length))
	v->data[loc] = gh_scm2double(val);
      else scm_misc_error(S_vct_setB,"invalid index",SCM_LIST3(obj,pos,val));
    }
  else scm_misc_error(S_vct_setB,"nil vct?",SCM_EOL);
  return(val);
}

static SCM vct_multiply(SCM obj1, SCM obj2)
{
  int i,lim;
  vct *v1,*v2;
  ERRVCT1(obj1,S_vct_multiplyB);
  ERRVCT2(obj2,S_vct_multiplyB);
  v1 = get_vct(obj1);
  v2 = get_vct(obj2);
  if ((v1) && (v2))
    {
      lim = MIN(v1->length,v2->length);
      for (i=0;i<lim;i++) v1->data[i] *= v2->data[i];
    }
  return(obj1);
}

static SCM vct_add(SCM obj1, SCM obj2)
{
  int i,lim;
  vct *v1,*v2;
  ERRVCT1(obj1,S_vct_addB);
  ERRVCT2(obj2,S_vct_addB);
  v1 = get_vct(obj1);
  v2 = get_vct(obj2);
  if ((v1) && (v2))
    {
      lim = MIN(v1->length,v2->length);
      for (i=0;i<lim;i++) v1->data[i] += v2->data[i];
    }
  return(obj1);
}

static SCM vct_scale(SCM obj1, SCM obj2)
{
  int i;
  vct *v1;
  float scl;
  ERRVCT1(obj1,S_vct_scaleB);
  ERRN2(obj2,S_vct_scaleB);
  v1 = get_vct(obj1);
  scl = gh_scm2double(obj2);
  if (v1)
    for (i=0;i<v1->length;i++) v1->data[i] *= scl;
  return(obj1);
}

static SCM vct_offset(SCM obj1, SCM obj2)
{
  int i;
  vct *v1;
  float scl;
  ERRVCT1(obj1,S_vct_offsetB);
  ERRN2(obj2,S_vct_offsetB);
  v1 = get_vct(obj1);
  scl = gh_scm2double(obj2);
  if (v1)
    for (i=0;i<v1->length;i++) v1->data[i] += scl;
  return(obj1);
}

static SCM vct_fill(SCM obj1, SCM obj2)
{
  int i;
  vct *v1;
  float scl;
  ERRVCT1(obj1,S_vct_fillB);
  ERRN2(obj2,S_vct_fillB);
  v1 = get_vct(obj1);
  scl = gh_scm2double(obj2);
  if (v1)
    for (i=0;i<v1->length;i++) v1->data[i] = scl;
  return(obj1);
}

static SCM list2vct(SCM lst)
{
  int len,i;
  vct *v;
  SCM scv;
  SCM_ASSERT(gh_list_p(lst),lst,SCM_ARG1,S_list2vct);
  len = gh_length(lst);
  scv = make_vct(len,(float *)CALLOC(len,sizeof(float)));
  v = get_vct(scv);
  for (i=0;i<len;i++) v->data[i] = (float)gh_scm2double(scm_list_ref(lst,gh_int2scm(i)));
  return(scv);
}

void init_vct(void)
{
#if HAVE_MAKE_SMOB_TYPE
  vct_tag = scm_make_smob_type_mfpe("vct",sizeof(vct),mark_vct,free_vct,print_vct,equalp_vct);
#else
  vct_tag = scm_newsmob(&vct_smobfuns);
#endif
  gh_new_procedure1_0(S_vct_length,vct_length);
  gh_new_procedure2_0(S_vct_ref,vct_ref);
  gh_new_procedure3_0(S_vct_setB,vct_set);
  gh_new_procedure2_0(S_vct_multiplyB,vct_multiply);
  gh_new_procedure2_0(S_vct_scaleB,vct_scale);
  gh_new_procedure2_0(S_vct_fillB,vct_fill);
  gh_new_procedure2_0(S_vct_addB,vct_add);
  gh_new_procedure2_0(S_vct_offsetB,vct_offset);
  gh_new_procedure1_0(S_make_vct,g_make_vct);
  gh_new_procedure1_0(S_vct_copy,copy_vct);
  gh_new_procedure1_0(S_vct_p,g_vct_p);
  gh_new_procedure1_0(S_list2vct,list2vct);
}
#endif