#include "ca.h"int (*ca_update_func)(ca_t *ca, int s, int e);/* get flat position in memory array from point   vector p */int ca_point(ca_t *ca, int *p) {  int i;  int run=1, c, n=0;  c = ca->cur;  for (i=0; i<(ca->dim); i++) {    n += run*p[i];    run *= ca->d[i];  }  return(n);}/* get point vector from flat memroy point p   and put into dest */int ca_vec(ca_t *ca, int *dest, int p) {  static int *d = NULL, *run = NULL, dim = 0;  int i, j, k;  if (dim!=(ca->dim)) {    if (d) free(d);    d = (int *)malloc(sizeof(int)*(ca->dim));    if (run) free(run);    run = (int *)malloc(sizeof(int)*(ca->dim));    dim = ca->dim;  }  run[0] = 1;  for (i=0; i<(dim-1); i++) {    run[i+1] = ca->d[i]*run[i];  }  for (i=(dim-1); i>=0; i--) {    dest[i] = p/run[i];    p -= (p/run[i])*run[i];  }}  int *ca_nei_addp(ca_t *ca, int *p, int *a) {  int i, j, k;  int carry=0;  for (i=0; i<(ca->dim); i++) {    p[i] += a[i] + ca->d[i];    p[i] %= ca->d[i];    /*    carry = ( (p[i]>ca->d[i]) ? 1 : 0 );    if (p[i]>=ca->d[i]) {      carry = p[i] - ca->d[i] + 1;    } else if (p[i]<0) {      carry = ca->d[i] + p[i];    } else carry=0;    p[i] += ca->d[i];    p[i] %= ca->d[i];    */  }  return(p);}/* get sum of neighbors for outer totalistic   CA */int ca_tot_sum(ca_t *ca, int p) {  int i, j, k=0;  int c;  int tfp;  static int *po=NULL, *t=NULL, dim=0;  if (dim!=ca->dim) {    dim = ca->dim;    if (po) {      free(po);      free(t);    }    po = (int *)malloc(sizeof(int)*(ca->dim));    t = (int *)malloc(sizeof(int)*(ca->dim));  }  c = ca->cur;  ca_vec(ca, t, p);  for (i=0; i<(ca->nei_num); i++) {    memcpy(po, t, sizeof(int)*dim);    ca_nei_addp(ca, po, ca->nei[i]);    tfp = ca_point(ca, po);    k += (ca->ca[c][tfp])*(ca->w[i]);    /* whoopsie, doesn't take into acount edges for > 1 dim */    /*    k += (ca->ca[c][(p + ca->fnei[i] + ca->fd)%ca->fd])*(ca->w[i]);    */  }  /*  if (ca->ca[c][p]) {    printf("k: %i (%i)\n", k, ca->wadd);  }  */  /* k is neighbor sum.  wadd is sum of weights.     since its outer totalistic, add the middle cell     value * sum of weights to get rule */  return(k + (ca->ca[c][p]*(ca->wadd)));}/* get next state for cell.   rule position = p_n p_(n-1) ... p_0 p */int ca_flat_state(ca_t *ca, int p) {  int i, j, k, c, n, *d;  int dig=0, digp=1, tfp;  static int *po=NULL, *t=NULL, dim=0;  if (dim!=ca->dim) {    dim = ca->dim;    if (po) {      free(po);      free(t);    }    po = (int *)malloc(sizeof(int)*(ca->dim));    t = (int *)malloc(sizeof(int)*(ca->dim));  }  c = ca->cur;  ca_vec(ca, t, p);  for (i=(ca->nei_num-1); i>=0; i--) {    memcpy(po, t, sizeof(int)*dim);    ca_nei_addp(ca, po, ca->nei[i]);    tfp = ca_point(ca, po);    dig += ca->ca[c][tfp]*digp;    digp *= ca->k;    /* whoops, doesn't handle edges in > 1 dimension */    /*    dig += ca->ca[c][ (ca->fnei[i] + p + ca->fd) % ca->fd ]*digp;    digp *= ca->k;    */  }  return(ca->r[dig]);}/* get rule position from flat position in memory */int ca_flat_rule(ca_t *ca, int p) {  int i, j, k, c, n, *d;  int dig, digp=1;  c = ca->cur;  dig = 0;  for (i=0; i<(ca->nei_num); i++) {    dig += ca->ca[c][ (ca->fnei[i] + p + ca->fd) % ca->fd ]*digp;    digp *= ca->dim;  }  /*  dig += ca->ca[c][p];  digp *= ca->dim;  for (i = ((ca->nei_num)/2); i<(ca->nei_num); i++) {    dig += ca->ca[c][ca->fnei[i]+p]*digp;    digp *= ca->dim;  }  */  return(dig);}int ca_tot_update(ca_t *ca, int s, int e) {  int i, j, c, n;  c = ca->cur;  n = (ca->cur + 1)%(ca->wlen);  for (i=s; i<=e; i++) {    ca->ca[n][i] = ca->tr[ca_tot_sum(ca, i)];  }  ca->cur++;  ca->cur%=ca->wlen;  return(0);}/*the way i'm doing it, there really isn't   any use of having s and e...*/int ca_flat_update(ca_t *ca, int s, int e) {  int i, j, c, n;  c = ca->cur;  n = (ca->cur+1)%(ca->wlen);  /* apply the rule to every cell */  for (i=s; i<=e; i++) {    ca->ca[n][i] = ca_flat_state(ca,i);  }  ca->cur++;  ca->cur%=ca->wlen;  return(0);}/* a more user friendly version of update.   ps start point, pe end point   DON'T USE (its not right)*/int ca_update(ca_t *ca, int *ps, int *pe) {  return(ca_flat_update(ca, ca_point(ca, ps), ca_point(ca, pe)));}/********************************* memory functions (alloc, init, free etc.)**********************************//* allocate the cell aut for the ca   assumed ca->dim & ca->d allocated and initialized */int ca_cell_alloc(ca_t *ca) {  int i, n=1;  /* allocate XxYxZ... size ca */  for (i=0; i<(ca->dim); i++) {    n*=ca->d[i];  }  //printf("ca: allocating ca (%i)\n", sizeof(int *)*(ca->wlen));  if (!(ca->ca = (int **)malloc(sizeof(int *)*(ca->wlen))))    return(-1);  for (i=0; i<(ca->wlen); i++) {    //printf("\t%i (%i)\n", i, sizeof(int)*n);    if (!(ca->ca[i] = (int *)malloc(sizeof(int)*n)))      return(-1);  }  ca->cur = 0;  ca->next = 1;  return(0);}ca_t *ca_alloc(int k, int dim, int wlen,	       int *d, int nei_num, int **nei,	       int *r, float **b, float *wb) {  ca_t *ca;  ca = (ca_t *)malloc(sizeof(ca_t));  return(ca_init(ca, k, dim, wlen, d, nei_num, nei, r, b, wb));}/* non totalistic init function */ca_t *ca_init(ca_t *ca, int k, int dim, int wlen,	      int *d, int nei_num, int **nei,	      int *r, float **b, float *wb) {  int i, j;  /*    ca_t *ca;    ca = (ca_t *)malloc(sizeof(ca_t));  */  ca->k = k;  ca->dim = dim;  ca->wlen = wlen;  ca->slen = wlen;  ca->go = 0;  /* dimensions of ca */  //printf("ca: allocating d (%i)\n", sizeof(int)*dim);  if (!(ca->d = (int *)malloc(sizeof(int)*dim)))    return(NULL);  ca->fd=1;  for (i=0; i<dim; i++) {    ca->d[i] = d[i];    ca->fd *= d[i];  }  ca->nei_num = nei_num;  /* neigbors */  //printf("ca: allocating nei & fnei (%i)\n", sizeof(int *)*nei_num);  if (!(ca->nei = (int **)malloc(sizeof(int *)*nei_num)))    return(NULL);  if (!(ca->fnei = (int *)malloc(sizeof(int *)*nei_num)))    return(NULL);  /* aloccate and assign relative neighbors */  for (i=0; i<nei_num; i++) {    //printf("\t%i (%i)\n", i, sizeof(int)*dim);    if (!(ca->nei[i] = (int *)malloc(sizeof(int)*dim)))      return(NULL);    for (j=0; j<dim; j++) {      ca->nei[i][j] = nei[i][j];    }  }  /* assign flat neighbors */  for (i=0; i<nei_num; i++) {    ca->fnei[i] = ca_point(ca, ca->nei[i]);  }  /* rule lookup */  ca->r_num = 1;  for (i=0; i<nei_num; i++) {    ca->r_num *= k;  }  //printf("ca: allocating r (%i)\n", sizeof(int)*(ca->r_num));  if (!(ca->r = (int *)malloc(sizeof(int)*(ca->r_num))))    return(NULL);  for (i=0; i<(ca->r_num); i++) {    ca->r[i] = r[i];  }  /* basis */  //printf("ca: allocating b (%i)\n", sizeof(float *)*dim);  if (!(ca->b = (float **)malloc(sizeof(float *)*dim)))    return(NULL);  for (i=0; i<dim; i++) {    //printf("ca:\t%i (%i)\n", i, sizeof(float)*4);    if (!(ca->b[i] = (float *)malloc(sizeof(float)*4)))      return(NULL);    for (j=0; j<dim; j++) {      ca->b[i][j] = b[i][j];    }    for (; j<4; j++) {      ca->b[i][j] = 0.0;    }  }  /* win vector */  //printf("ca: allocating wb (%i)\n", sizeof(float)*(dim+1));  if (!(ca->wb = (float *)malloc(sizeof(float)*(dim+1))))    return(NULL);  for (i=0; i<(dim+1); i++) {    ca->wb[i] = wb[i];  }  ca->delay = 1000;  ca->cur = 0;  ca->next = 1;  if (ca_cell_alloc(ca)<0) {    return(NULL);  }  for (j=0; j<ca->wlen; j++) {    for (i=0; i<ca->fd; i++) {      ca->ca[j][i] = 0;    }  }  ca->diff = 0;  ca->slen = 1;  return(ca);}/* outer totalistic init */ca_t *ca_tot_init(ca_t *ca, int k, int dim, int wlen,		  int *d, int nei_num, int **nei, int *w,		  int *r, float **b, float *wb) {  int i, j;  /*    ca_t *ca;    ca = (ca_t *)malloc(sizeof(ca_t));  */  ca->k = k;  ca->dim = dim;  ca->wlen = wlen;  ca->slen = wlen;  ca->go = 0;  /* dimensions of ca */  if (!(ca->d = (int *)malloc(sizeof(int)*dim)))    return(NULL);  ca->fd=1;  for (i=0; i<dim; i++) {    ca->d[i] = d[i];    ca->fd *= d[i];  }  ca->nei_num = nei_num;  /* neigbors */  if (!(ca->nei = (int **)malloc(sizeof(int *)*nei_num)))    return(NULL);  if (!(ca->fnei = (int *)malloc(sizeof(int *)*nei_num)))    return(NULL);  /* aloccate and assign relative neighbors */  for (i=0; i<nei_num; i++) {    if (!(ca->nei[i] = (int *)malloc(sizeof(int)*dim)))      return(NULL);    for (j=0; j<dim; j++) {      ca->nei[i][j] = nei[i][j];    }  }  /* weights */  ca->w = (int *)malloc(sizeof(int)*nei_num);  for (i=0; i<nei_num; i++) {    ca->w[i] = w[i];  }  ca->wadd=1;  for (i=0; i<ca->nei_num; i++) {    ca->wadd += ca->w[i];  }  /* assign flat neighbors */  for (i=0; i<nei_num; i++) {    ca->fnei[i] = ca_point(ca, ca->nei[i]);  }  /* outer totalistic rule lookup */  ca->tr_num = (nei_num+1)*k;  if (!(ca->tr = (int *)malloc(sizeof(int)*(ca->tr_num))))    return(NULL);  for (i=0; i<(ca->tr_num); i++) {    ca->tr[i] = r[i];  }  /* basis */  if (!(ca->b = (float **)malloc(sizeof(float *)*dim)))    return(NULL);  for (i=0; i<dim; i++) {    if (!(ca->b[i] = (float *)malloc(sizeof(float)*dim)))      return(NULL);    for (j=0; j<dim; j++) {      ca->b[i][j] = b[i][j];    }  }  /* win vector */  if (!(ca->wb = (float *)malloc(sizeof(float)*(dim+1))))    return(NULL);  for (i=0; i<(dim+1); i++) {    ca->wb[i] = wb[i];  }  ca->delay = 1000;  ca->cur = 0;  ca->next = 1;  if (ca_cell_alloc(ca)<0) {    return(NULL);  }  ca->slen = 1;  ca->diff = 0;  return(ca);}void ca_free(ca_t *ca) {  int i, j, k;  //printf("ca: freeing d\n");  free(ca->d);  for (i=0; i<ca->nei_num; i++) {    //printf("ca:\t%i\n", i);    free(ca->nei[i]);  }  //printf("ca: freeing nei & fnei\n");  free(ca->nei);  free(ca->fnei);  for (i=0; i<ca->wlen; i++) {    //printf("ca:\t%i\n", i);    free(ca->ca[i]);  }  //printf("ca: freeing ca\n");  free(ca->ca);  if (ca->w) {    //printf("ca: freeing w\n");    free(ca->w);  }  for (i=0; i<ca->dim; i++) {    //printf("ca:\t%i\n", i);    free(ca->b[i]);  }  //printf("ca: freeing b\n");  free(ca->b);  if (ca->wb) {    //printf("ca: freeing wb\n");    free(ca->wb);  }  if (ca->r) {    //printf("ca: freeing r\n");    free(ca->r);  }  if (ca->tr) {    //printf("ca: freeing tr\n");    free(ca->tr);  }  //printf("ca: freeing mother ca\n");  free(ca);}/**************************************