#include <stdio.h>
#include <math.h>
#include <mpl.h>
#include <mpml.h>
#include "Wempar.h"

static int first_time_cos=0;
static int first_time_sin=0;
static int pre_N=0;
static WmdMat *COS_X;
static WmdMat *SIN_X;

static void
  Wmp_fft_init(N,gcos,odd)
int N,gcos,odd;
{
  int nA,nB,irow,icol,TWO_N;
  double value;
  WmdVex *cos_array;

  nA = N >> 1;
  if (odd) nB = nA;
  else     nB = nA + 1;
  TWO_N=2*N;
 
  if (N == 0) { /* free-up memmory */
    WmdMat_free(COS_X);
    WmdMat_free(SIN_X);
    pre_N=0;
    return;
  }

  if (pre_N != N) {
    first_time_sin=first_time_cos=0;
    pre_N=N;
  }

  if ((first_time_cos == 0) && (gcos == 1)) { /* if not initialize, do so */
    value=M_PI/N;
    cos_array=(WmdVex *)WmdVex_open(TWO_N);
    WmdVex_fill_linear(cos_array,value);
    WmdVex_func(cos_array,cos_array,p_cos);
    COS_X=(WmdMat *)WmdMat_open(nB,nA);
    for (irow=0; irow < COS_X->hdr.nrow; irow++) {
      for (icol=0; icol < COS_X->hdr.ncol; icol++) {
	value=WmdVex_get(cos_array,((2*irow+odd)*icol)%TWO_N);
	WmdMat_put(COS_X,irow,icol,value);
      }
    }
    first_time_cos=1;
    WmdVex_free(cos_array);
  }

  if ((first_time_sin == 0) && (gcos == 0)) { /* if not initialize, do so */
    value=M_PI/N;
    cos_array=(WmdVex *)WmdVex_open(TWO_N);
    WmdVex_fill_linear(cos_array,value);
    WmdVex_func(cos_array,cos_array,p_sin);
    SIN_X=(WmdMat *)WmdMat_open(nB,nA);
    for (irow=0; irow < SIN_X->hdr.nrow; irow++) {
      for (icol=0; icol < SIN_X->hdr.ncol; icol++) {
	value=WmdVex_get(cos_array,((2*irow+odd)*icol)%TWO_N);
	WmdMat_put(SIN_X,irow,icol,value);
      }
    }
    first_time_sin=1;
    WmdVex_free(cos_array);
  }
}

void
  Wmp_fft(A,B,gcos,odd,N)
WmdVex *A;
WmdVey *B;
int gcos,odd,N;
{
  /* initalize DFT cofficients */
  Wmp_fft_init(N,gcos,odd);
  if (N <  1) return;

  if (gcos) /* do a cosine seires only */
    WmdMat_vexmul(COS_X,A,B);
  else /* do a sine series only */
    WmdMat_vexmul(SIN_X,A,B);
}

static double *cos_arr;

#define PI M_PI

void
  WmFFT_init(gcos,N)
int gcos,N;
{
  register double rtmp;
  plural double dcos;
  int i,N4;

  if (N < 1) return;

  N4 = 2*N;
  rtmp = PI/N;
  if (cos_arr!=NULL) { free( (char *)cos_arr); cos_arr=NULL;}
  cos_arr = (double *)mp_new(N4*sizeof(double)); 
  if (N4 > nproc) {
    if (gcos) for (i=0; i < N4; i++)   *(cos_arr+i) = cos(i*rtmp);
    else      for (i=0; i < N4; i++)   *(cos_arr+i) = sin(i*rtmp);
  }
  else {
    if (gcos) {if (iproc < N4)     dcos=p_cos(iproc*rtmp);}
    else      {if (iproc < N4)     dcos=p_sin(iproc*rtmp);}
    for (i=0; i < N4; i++)       *(cos_arr+i)=proc[i].dcos;
  }
}

void
  WmFFT_free()
{
  if (cos_arr) { free(cos_arr);  cos_arr=NULL; }
}

void
  WmFFT_do(A,gcos,odd,sym,N,n)  
plural double *A; 
int gcos,odd,sym,N,n;
{ 
  if (N==0) {
    WmFFT_free();
    return;
  }
  
  WmFFT_init(gcos,N);

  if (sym) cs_dft_sym (A,odd,N,n);
  else     cs_dft_full(A,odd,N,n);

  WmFFT_free();
} /*cs_fft*/

/*
 * Do a full cos/sin series (0..N/2), non-symmetric case.                       
 *  Ac or As pts to N/2 items and returns N/2+1 (!odd) or N/2 (odd) items.
 *
 */
cs_dft_full(A,odd,N0,np)  
plural double A[];
int odd,N0,np;
{ 
  plural double *B;
  int N2,N3,N4;
  plural int i,k,n;
  double rtmp;
  register plural double rrtmp;

  odd = !!odd;
  N2 = N0>>1;  
  N3 = N2 + !(odd);  
  N4 = 2*N0;

  if (np > nproc) Wmp_bugs("cs_dft_full","Not enough processors avaliable");
  /*Check that enough memory is allocated for B*/
  B = (plural double *)p_new( N2*sizeof(double) ); 
  if (iproc < np) 
    for (n=0; n<N2; n++) *(B+n) = *(A+n);

  /*Do the dft*/
  if (iproc < np) {
    for (i=0; i<N3; i++) {
      rrtmp = 0.0;  
      k =  2*i + odd;
      for (n=0; n<N2; n++) 
	rrtmp += B[n] * cos_arr[(n*k)%N4]; 
      *(A+i) = rrtmp;
    }
  }
  p_free(B);

} /*cs_dft_full*/

/*
 * Do a full cos/sin series (0..N/2), symmetric case.                       
 *  Ac or As pts to N/2 items and returns N/2+1 (!odd) or N/2 (odd) items.
 *
 */
cs_dft_sym(A,odd,N0,np)
plural double *A;  
int odd,N0,np;
{
  plural double *B;
  plural int i,k,n;
  int N2,N3,N4,NN4;
  register plural double rrtmp;

  odd = !!odd;
  N2 = N0>>1; 
  N3 = ((!odd) && (!(N2&1))) + ((N2+1)>>1);
  N4 = 2*N0; 
  NN4 = N2>>1;

  /*Check that enough memory is allocated for B*/
  B = (plural double *)p_new( NN4*sizeof(double) ); 
  if (iproc < np) 
    for (n=0; n<NN4; n++) *(B+n) = *(A+n);

  /*Do the dft*/
  if (iproc < np) {
    for (i=0; i<N3; i++) {
      rrtmp = 0.0;  
      k =  2*i + odd;
      for (n=0; n<N2; n++) 
	if (n&1) rrtmp += *(B+(n>>1)) * *(cos_arr + (n*k)%N4); 
      *(A+i) = rrtmp;
    }
  }
  p_free(B);
} /*cs_dft_sym*/






