/*======================================================*/
/*
K-K transformation program 
HP version         K.Kobayshi 1996. 8. 26  
*/
/*======================================================*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h> /** important!! **/
#include </home/nakao/c/include/n_string.h>

#define N 32768

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(data,nn,isign)
     float data[];
     int nn,isign;
{
  int n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  float tempr,tempi;
  
  n=nn << 1;
  j=1;
  for (i=1;i<n;i+=2) {
    if (j > i) {
      SWAP(data[j],data[i]);
      SWAP(data[j+1],data[i+1]);
    }
    m=n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax=2;
  while (n > mmax) {
    istep=2*mmax;
    theta=6.28318530717959/(isign*mmax);
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for (m=1;m<mmax;m+=2) {
      for (i=m;i<=n;i+=istep) {
	j=i+mmax;
	tempr=wr*data[j]-wi*data[j+1];
	tempi=wr*data[j+1]+wi*data[j];
	data[j]=data[i]-tempr;
	data[j+1]=data[i+1]-tempi;
	data[i] += tempr;
	data[i+1] += tempi;
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
  }
}

#undef SWAP

void realft(data,n,isign)
     float data[];
     int n,isign;
{
  int i,i1,i2,i3,i4,n2p3;
  float c1=0.5,c2,h1r,h1i,h2r,h2i;
  double wr,wi,wpr,wpi,wtemp,theta;
  void four1();
  
  theta=3.141592653589793/(double) n;
  if (isign == 1) {
    c2 = -0.5;
    four1(data,n,1);
  } else {
    c2=0.5;
    theta = -theta;
  }
  wtemp=sin(0.5*theta);
  wpr = -2.0*wtemp*wtemp;
  wpi=sin(theta);
  wr=1.0+wpr;
  wi=wpi;
  n2p3=2*n+3;
  for (i=2;i<=n/2;i++) {
    i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
    h1r=c1*(data[i1]+data[i3]);
    h1i=c1*(data[i2]-data[i4]);
    h2r = -c2*(data[i2]+data[i4]);
    h2i=c2*(data[i1]-data[i3]);
    data[i1]=h1r+wr*h2r-wi*h2i;
    data[i2]=h1i+wr*h2i+wi*h2r;
    data[i3]=h1r-wr*h2r+wi*h2i;
    data[i4] = -h1i+wr*h2i+wi*h2r;
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
  }
  if (isign == 1) {
    data[1] = (h1r=data[1])+data[2];
    data[2] = h1r-data[2];
  } else {
    data[1]=c1*((h1r=data[1])+data[2]);
    data[2]=c1*(h1r-data[2]);
    four1(data,n,-1);
  }
}

void sinft(y,n)
     float y[];
     int n;
{
  int j,m=n/2,n2=n+2;
  float sum,y1,y2;
  double theta,wi=0.0,wr=1.0,wpi,wpr,wtemp;
  void realft();
  
  theta=3.14159265358979/(double) n;
  wtemp=sin(0.5*theta);
  wpr = -2.0*wtemp*wtemp;
  wpi=sin(theta);
  y[1]=0.0;
  for (j=2;j<=m+1;j++) {
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
    y1=wi*(y[j]+y[n2-j]);
    y2=0.5*(y[j]-y[n2-j]);
    y[j]=y1+y2;
    y[n2-j]=y1-y2;
  }
  realft(y,m,1);
  y[1]*=0.5;
  sum=y[2]=0.0;
  for (j=1;j<=n-1;j+=2) {
    sum += y[j];
    y[j]=y[j+1];
    y[j+1]=sum;
  }
}

void cosft(y,n,isign)
     float y[];
     int n,isign;
{
  int j,m,n2;
  float enf0,even,odd,sum,sume,sumo,y1,y2;
  double theta,wi=0.0,wr=1.0,wpi,wpr,wtemp;
  void realft();
  
  theta=3.14159265358979/(double) n;
  wtemp=sin(0.5*theta);
  wpr = -2.0*wtemp*wtemp;
  wpi=sin(theta);
  sum=y[1];
  m=n >> 1;
  n2=n+2;
  for (j=2;j<=m;j++) {
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
    y1=0.5*(y[j]+y[n2-j]);
    y2=(y[j]-y[n2-j]);
    y[j]=y1-wi*y2;
    y[n2-j]=y1+wi*y2;
    sum += wr*y2;
  }
  realft(y,m,1);
  y[2]=sum;
  for (j=4;j<=n;j+=2) {
    sum += y[j];
    y[j]=sum;
  }
  if (isign == -1) {
    even=y[1];
    odd=y[2];
    for (j=3;j<=n-1;j+=2) {
      even += y[j];
      odd += y[j+1];
    }
    enf0=2.0*(even-odd);
    sumo=y[1]-enf0;
    sume=(2.0*odd/n)-sumo;
    y[1]=0.5*enf0;
    y[2] -= sume;
    for (j=3;j<=n-1;j+=2) {
      y[j] -= sumo;
      y[j+1] -= sume;
    }
  }
}

/*imtre.c*/
main(argc,argv)
int argc;
char **argv;
{
  float a[N],b[N],x,y;
  int i,j,m;
  FILE *ifp,*ofp,*fopen();
  char fname[25];

  if(argc!=3){
    printf("\n%s [infile] [outfile] \n\n",argv[0]);
    exit(0);
  }
  
  if((ifp=fopen(argv[1],"r"))==NULL){
    printf("This file (%s) can not open !!\n",argv[1]);
    exit(0);
  }
  if((ofp=fopen(argv[2],"w"))==NULL){
    printf("This file (%s) can not open !!\n",argv[2]);
    exit(0);
  }
  
  printf("Enter data number ");
  scanf("%d",&m);
  
  for ( i=1; i <= m; ++i) { 
    fscanf(ifp,"%g %g", &a[i], &b[i]); 
  }
  
  for(i=m+1;i<=N;++i){
    a[i]=0.0;
    b[i]=0.0;
  }

  fclose(ifp);
  sinft(b,N);
  
  for (i=1; i <= N; ++i) {
    b[i]*=2.0/N;
  } 

  cosft(b,N,1);

  for ( i=1; i<=m; ++i) {
    fprintf(ofp, "%g %g\n",a[i],-1*b[i]);
  }

  printf("finished!!\n");
  fclose(ofp);
}