#ifndef FAST_H
#define FAST_H

#define PROD5(A,B,R) \
{ typeof (R) __x,__y,__t;	\
  typeof (&A[0]) __pa=A;	\
  typeof (&B[0]) __pb=B;	\
                          __x=*__pa++;__y=*__pb++;	\
  R=R-R;    __t=__x*__y;  __x=*__pa++;__y=*__pb++;	\
  R+=__t;   __t=__x*__y;  __x=*__pa++;__y=*__pb++;	\
  R+=__t;   __t=__x*__y;  __x=*__pa++;__y=*__pb++;	\
  R+=__t;   __t=__x*__y;  __x=*__pa++;__y=*__pb++;	\
  R+=__t;   __t=__x*__y;	\
  R+=__t;}


/* Autocorrelation : R[0:K] is autocorrelation of X[M:N-1]  i.e.
      R[k] = Sum X[i]*X[i-k]   for    M<=i<N     */

#define AUTOCORR(X,R,K,M,N) \
{\
    int  ii,jj;\
    for(ii=0; ii<=K; ii++) {\
	R[ii] = 0;\
	for(jj=M; jj<N; jj++) R[ii] += X[jj]*X[jj-ii];}}

/* if AC is autocorrelation of X (as in above) , then
   R[i] = D * R[i] + AC[i]    */

#define RECACORR(X,R,K,M,N,D)	\
{   int __k,__i,__M1=M+1;	\
    typeof (&R[0]) __rp=R;	\
    for(__k=0; __k<=K; __k++) {	\
	typeof (*R) __corr, __t;	\
	typeof (&X[0]) __xip=X+M, __xkp=__xip-__k;	\
	typeof(X[0]) __xi,__xk;	\
	__corr=__corr-__corr;	\
	__xi=*__xip++;	\
	__xk=*__xkp++;	\
	for(__i=__M1;__i<N;__i++) {	\
	    __t=__xi*__xk; __xi=*__xip++;	\
	    __corr+=__t; __xk=*__xkp++;	\
	}	\
	__t=__xi*__xk;	\
	__corr+=__t;	\
	__corr=*__rp*(D)+__corr;	\
	*__rp++=__corr;}}


/** Z[i] = X[i]*Y[i] **/

/** Hacked this so that it only works with reals. This is necessary
    because the non-GNU C compiler has trouble with typeof() **/
#define VPROD(X,Y,Z,N)  \
{   real* __xp=X;       \
    real* __yp=Y;      \
    real* __zp=Z;      \
    int i;      \
    real __xi=*__xp++; \
    real __yi=*__yp++; \
    real __zi; \
    for(i=1;i<N;i++) {  \
        __zi=__xi*__yi; \
        __xi=*__xp++;   \
        __yi=*__yp++;   \
        *__zp++=__zi;}  \
    __zi=__xi*__yi;     \
    *__zp=__zi;}


/** Hacked this so that it only works with reals */
#define RSHIFT(START,LENGTH,SHIFT) \
{ real* __from = &(START)[LENGTH-1]; \
  real* __to = &(START)[LENGTH+SHIFT-1]; \
  int __i; for(__i=0; __i<LENGTH; __i++)      *__to-- = *__from--;\
}


/** Must have at least two elements, i.e. N >= 2.
  computes sum(x[i]*y[i]) **/

/** Hacked this so that it only works with reals. This is necessary
    because the non-GNU C compiler has trouble with typeof() **/
#define DOTPROD(X,Y,N,Z) \
{ int __i;	\
  real __r=0, __t;	\
  real __x1;	\
  real __y1;	\
  real* __xp = &(X)[0]; \
  real* __yp = &(Y)[0]; \
                              __x1=*__xp++; __y1=*__yp++; 	\
             __t=__x1*__y1;   __x1=*__xp++; __y1=*__yp++; 	\
   for(__i=0; __i<(N)-2; __i++) {	\
   __r+=__t; __t=__x1*__y1;   __x1=*__xp++; __y1=*__yp++; }	\
   __r+=__t; __t=__x1*__y1;   	\
   __r+=__t; *Z=__r;}

#endif /*FAST_H*/