#define NPOLE 6
#define NPOLEM1 5
float cq,outold;
main()
{
	float z,sig[100],past[100],d,c[100],out[100],x;
	float bwarppol(),past4[100],out2[100],ballpole(),allpole();
	float warpset(),cl,warppol(),zz,past2[100],past3[100];
	int i,npoles,nvals,j,jcount;
	jcount = 0;
	for(j=0; j<100; j++) {
	sig[j] = past[j] = out[j]= past2[j] = past3[j] = past4[j] = out2[j];
	c[j] = 1./(float)((j+1)*100);
	}
	z=d=0;
	warpinit();
	cl=warpset(d,c);
	x=1;
	sig[0]=1;
	ballpole(sig,&jcount,NPOLE,past3,c,out,30);
	bwarppol(sig,NPOLE,past4,d,c,out2,30);
	warpinit();
	cl=warpset(d,c);
	jcount=0;
	for(i=0; i<30; i++) {
		z=allpole(x,&jcount,NPOLE,past,c);
		zz=warppol(x,past2,d,c);
		printf("%d %e %e %e %e %e\n",i,z,zz,cl,out[i],out2[i]);
		x=0;
	}
}
float allpole(x,jcount,npoles,past,c)
float x,*past,*c;
long  *jcount,npoles;
{
	int j,nfint;
	for(j= *jcount, nfint=0;  nfint<npoles;  nfint++,j++)
		x += (*(c+nfint) * *(past+j));
	*(past+ *jcount) = *(past+*jcount+npoles) = x;
	*jcount = (*jcount + 1) % npoles;
	return(x);
}
float warppol(sig,past,d,c)
float sig,*past,d,*c;
{
	float temp1,temp2;
	register n;

	temp1 = past[NPOLEM1];
	past[NPOLEM1] = cq * outold - d * past[NPOLEM1];
	for(n=NPOLE-2; n>=0; n--) {
		temp2 = past[n];
		past[n] = d * (past[n+1] - past[n]) + temp1;
		temp1 = temp2;
		}
	for(n=0;n<NPOLE;n++)  sig += c[n] * past[n];
	outold = sig;
	return(sig);
}

float warpset(d,c)     
float d,*c;
{
	register m;
	float cl;

	for(m=1; m<NPOLE; m++)   c[m] += d * c[m-1];
	cl = 1./(1.-d * c[NPOLEM1]);
	cq = cl * (1. - d * d);
	return(cl);
}

warpinit()
{
	outold = 0;
}


float ballpole(x,jcount,npoles,past,c,out,nvals)
float *x,*past,*c,*out;
long  *jcount,npoles,nvals;
{
	register i,j,nfint;
	register float *temp;
	for(i=0;i<nvals;++i){
		*temp = *x++;
		for(j= *jcount, nfint=0;  nfint<npoles;  nfint++,j++)
			*temp += (*(c+nfint) * *(past+j));
		*out++ =  *(past+ *jcount) = *(past+*jcount+npoles) = *temp;
		*jcount = (*jcount + 1) % npoles;
		}
}
float bwarppol(sig,npoles,past,d,c,out,nvals)
float *sig,*past,d,*c,*out;
{
	float temp1,temp2;
	int i,n;
	for(i=0; i<nvals; i++) {
		temp1 = past[npoles-1];
		past[npoles-1] = cq * outold - d * past[npoles-1];
		for(n=npoles-2; n>=0; n--) {
			temp2 = past[n];
			past[n] = d * (past[n+1] - past[n]) + temp1;
			temp1 = temp2;
			}
		for(n=0;n<npoles;n++)   *sig += c[n] * past[n];
		*out++ = outold = *sig++;
	}
}