/* nag_fft_real(c06eac) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagc06.h>

#define NMAX 20

int main(void)
{
  Integer j, n, n2, nj;
  double a[NMAX], b[NMAX], x[NMAX], xx[NMAX];

  Vprintf("c06eac Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  while (scanf("%ld", &n)!=EOF)
    if (n>1 && n<=NMAX)
      {
        for (j = 0; j<n; j++)
          {
            Vscanf("%lf", &x[j]);
            xx[j] = x[j];
          }
        /* Calculate transform */
        c06eac(n, x, NAGERR_DEFAULT);
        /* Calculate full complex form of Hermitian result */
        a[0] = x[0];
        b[0] = 0.0;
        n2 = (n-1)/2;
        for (j = 1; j<=n2; j++)
          {
            nj = n - j;
            a[j] = x[j];
            a[nj] = x[j];
            b[j] = x[nj];
            b[nj] = -x[nj];
          }
        if (n % 2==0)
          {
            a[n2+1] = x[n2+1];
            b[n2+1] = 0.0;
          }
        Vprintf("\nComponents of discrete Fourier transform\n");
        Vprintf("\n         Real       Imag \n\n");
        for (j = 0; j<n; j++)
          Vprintf("%3ld %10.5f %10.5f\n", j, a[j], b[j]);
        /* Calculate inverse transform */
        c06gbc(n, x, NAGERR_DEFAULT);
        c06ebc(n, x, NAGERR_DEFAULT);
        Vprintf("\nOriginal sequence as restored by inverse transform\n");
        Vprintf("\n       Original   Restored\n\n");
        for (j = 0; j<n; j++)
          Vprintf("%3ld %10.5f %10.5f\n", j, xx[j], x[j]);
      }
    else
      {
        Vfprintf(stderr,"Invalid value of n\n");
        return EXIT_FAILURE;
      }
  return EXIT_SUCCESS;
}