#include <stdio.h>
#include <math.h>
#define GETEXP(a,b) var.d = 2*a; var.i >>= 23; b = 127-(var.i & 127)
#define MAX(a,b) ((a)>=(b)) ? (a) : (b)
#define MM(m,i) mm[i+m*128]

/* globals  */
union Var_
{
   float d;
   int i;
};
union Var_ var;

/*=================================================================
  Computes the correlation dimension, after P. Grassberger.
 *================================================================*/

void correl(x, kmax, eps, logimax, membed, mindel, mm)
float *x;
float eps;
int kmax, membed, mindel, logimax;
int *mm;
{
  int box[256][256];
  int llist[100000];
  int imax1, m, i, j, k;
  int i1, j1, i2, j2, l1, l2;
  int kp, kpn, logdx;
  float epsinv, x1, x2;
  float dx,ddx, dx2;

  imax1 = 1;
  imax1 <<= logimax;
  imax1 -= 1;
  epsinv=(float) 1./eps;

  for (m=0;m<=membed;m++)
    for(i=0;i<128;i++)
      MM(m,i)=0;

  for(i=0;i<imax1;i++)
    for(j=0;j<imax1;j++)
      box[i][j]=0;

  for(i=0;i<=kmax;i++)
    llist[i]=0;

  x1=x[0];
  i1 = (int)(x1*epsinv);
  i1 = i1 & imax1;

  for(k=1;k<=kmax-membed+1;k++)
    {
      x2=x[k];
      i2=(int) (x2*epsinv);
      i2=i2 & imax1;
      for(j1=i1-1;j1<=i1+1;j1++)
	{
	  l1=j1&imax1;
	  for(j2=i2-1;j2<=i2+1;j2++)
	    {
	      l2=j2&imax1;
	      kp=box[l1][l2];
	      while((kp>0)&&(kp<(k-mindel)))
		{ 
		  dx=fabs(x1-x[kp-1]);
		  if(dx<eps)
		    {
		      dx2=fabs(x2-x[kp]);
		      m=2;
		      while(dx2<eps && m <= membed)
			{
			  dx=MAX(dx,dx2);
			  if(dx==0)
			    logdx=127;
			  else
			      GETEXP(dx,logdx);
/*			  fprintf(stderr,"%f %d\n",dx,logdx); */
			  MM(m,logdx)++;
/*			  fprintf(stderr,"%d\t%d\t%d\n",m,logdx,MM(m,logdx)); */
			  if(m<membed)
			    {
			      m++;
			      dx2=fabs(x[k+m-2]-x[kp+m-2]);
			    }
			  else
			    m++;
			}

		    }
		  kp=llist[kp];
		}
	    }
	}
      kpn=box[i1][i2];
      if(kpn==0)
	box[i1][i2]=k;
      else
	{
	  while(kpn>0)
	    {
	      kp=kpn;
	      kpn=llist[kp];
	    }
	  llist[kp]=k;
	}
      x1=x2 ;
      i1=i2;
    }
  
  for(m=2;m<=membed;m++)
    for(i=126;i>=0;i--)
      {
	MM(m,i)+=MM(m,i+1);
/*	fprintf(stderr,"%d\t%d\t%d\n",m,i,MM(m,i)); */
      }
}