/*
  Combined Periodic Filter Transform on $2q$ Samples

  pcqfilter( out[], in[], q, h[], g[], L):
  [0] For n=0 to q-1, do [1] through [4]
  [1]    Let out[n] = 0, let out[n+q] = 0
  [2]    For k=0 to L-1, do [3] and [4]
  [3]       Accumulate  out[n]  += h[k]*in[(k+2*n)%(2*q)]
  [4]       Accumulate out[n+q] += g[k]*in[(k+2*n)%(2*q)]
  [5] Return


  Inverse of the Combined Periodic Filter Transform on $2q$ Samples

  ipcqfilter( out[], in[], q, h[], g[], L):
  [0] For k2=0 to q-1, do [1] to [4+]
  [1]    Let out[2*k2] = 0, let out[2*k2+1] = 0
  [2]    For n2=0 to L/2-1, do [3] to [4+]
  [3]       Accumulate  out[2*k2]  +=  h[2*n2]*in[(k2-n2) mod q]
  [3+]      Accumulate  out[2*k2]  +=  g[2*n2]*in[((k2-n2) mod q) + q]
  [4]       Accumulate out[2*k2+1] += h[2*n2+1]*in[(k2-n2) mod q]
  [4+]      Accumulate out[2*k2+1] += g[2*n2+1]*in[((k2-n2) mod q) + q]
  [5] Return

*/

/* Compute x%M for integer M>1 and any int x */
int
mod(int x, int M)
{
  if(x<0) x-= x*M;	/* x-x*modulus>0 is congruent to x mod M */
  return x%M;
}

void
pcqfilter(float out[], float in[], int q,
	  const float h[], const float g[], int L)
{
  int n, k;
  for(n=0; n<q; n++) {
    out[n]=out[n+q]=0;
    for(k=0;k<L;k++) {
      out[n] += h[k]*in[(k+2*n)%(2*q)];
      out[n+q] += g[k]*in[(k+2*n)%(2*q)];
    }
  }
  return;
}

void
ipcqfilter(float out[], float in[], int q,
	   const float h[], const float g[], int L)
{
  int n2, k2;
  for(k2=0; k2<q; k2++) {
    out[2*k2]=out[2*k2+1]=0;
    for(n2=0;n2<L/2;n2++) {
      out[2*k2] += h[2*n2]*in[mod(k2-n2, q)];
      out[2*k2] += g[2*n2]*in[mod(k2-n2, q) + q];
      out[2*k2+1] += h[2*n2+1]*in[mod(k2-n2, q)];
      out[2*k2+1] += g[2*n2+1]*in[mod(k2-n2, q) + q];
    }
  }
  return;
}