#define _USE_MATH_DEFINES
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#ifndef M_PI
#define M_PI 3.141592653589793
#endif

typedef double metric_wave_f(double);
typedef double metric_distance_f(double,double); // metric
typedef double metric_phase_f(double,double);

#define METRICS \
  X(euclidian) \
  X(taxicab)

enum metric_e {
#define X(X) WF_ ## X,
METRICS
#undef X
};

#define X(X) \
  metric_wave_f     metric_wave_ ## X; \
  metric_distance_f metric_distance_ ## X; \
  metric_phase_f    metric_phase_ ## X;
METRICS
#undef X

typedef struct metric {
  const char* name;
  metric_wave_f*     f_wave; // This should return 1 for input 0, and 0 for 0.25
  metric_distance_f* f_distance;
  metric_phase_f*    f_phase;
} metric_t;

const metric_t metric_list[] = {
#define X(X) { \
    .name=#X, \
    .f_wave     = metric_wave_ ## X, \
    .f_distance = metric_distance_ ## X, \
    .f_phase    = metric_phase_ ## X, \
  },
  METRICS
#undef X
};
enum { metric_count = sizeof(metric_list)/sizeof(*metric_list) };

double periodic(double t){
  t = fmod(t,1);
  if(t < 0) t += 1;
  return fabs(t);
}

// euclidian metric
double metric_wave_euclidian(double t){
  return cos(2.*M_PI*t);
}
double metric_distance_euclidian(double x, double y){
  return sqrt(x*x + y*y);
}
double metric_phase_euclidian(double x, double y){
  return atan2(y,x) / (2*M_PI);
}

// taxicab metric
double metric_wave_taxicab(double t){
  t = periodic(t);
  return fabs(t*2-1)*2-1;
}
double metric_distance_taxicab(double x, double y){
  return fabs(x) + fabs(y);
}
double metric_phase_taxicab(double x, double y){
  const bool xp = x<0;
  const bool yp = y<0;
  const bool bx = xp != yp;
  const double distance = metric_distance_taxicab(x,y);
  int quadrant = 2*yp + bx;
  if(bx){
    return (quadrant + fabs(x)/distance) / 4.;
  }else{
    return (quadrant + fabs(y)/distance) / 4.;
  }
}


void discrete_fourier_transform_sc(const metric_t*const metric, size_t n, double result[static restrict n][2], const double input[static n]){
  memset(result, 0, sizeof(double[n][2]));
  for(size_t i=0; i<n; i++){
    double v = input[i];
    for(size_t f=0; f<n; f++){
      double t = (double)(i*f)/n;
      double f1 = metric->f_wave(t);
      double f2 = metric->f_wave(t+0.25);
      double c1 = v * f1 / n;
      double c2 = v * f2 / n;
      result[f][0] += c1;
      result[f][1] += c2;
    }
  }
}

#define EPSILON 0.0000001

void discrete_fourier_transform(const metric_t*const metric, size_t n, double result[static n][2], const double input[static restrict n]){
  discrete_fourier_transform_sc(metric, n, result, input);
  for(size_t i=0; i<n; i++){
    double frequency = metric->f_distance(result[i][0], result[i][1]);
    double phase = metric->f_phase(result[i][0], result[i][1]);
    if(frequency < EPSILON && frequency > -EPSILON){
      frequency = 0;
      phase = 0;
    }
    if(phase+EPSILON <  0) phase += 1;
    result[i][0] = frequency;
    result[i][1] = phase; // 0..1
  }
}

void inverse_discrete_fourier_transform(const metric_t*const metric, size_t n, double result[static n], const double input[static restrict n][2]){
  memset(result, 0, sizeof(double[n]));
  for(size_t o=0; o<n; o++){
    double v = 0;
    for(size_t f=0; f<n; f++){
      double t = (double)(o*f)/n + input[f][1];
      v += metric->f_wave(t) * input[f][0];
    }
    result[o] = v;
  }
}

void pnum(double x){
  char b[42];
  int n=snprintf(b,42," % 10.6lf", x);
  while(--n>0 && b[n] == '0')
    b[n] = ' ';
  if(n > 0 && b[n] == '.')
    b[n] = ' ';
  if(!strcmp(b,"  -0       "))
    b[2] = ' ';
  printf("%s", b);
}

int main(int argc, char* argv[]){
  const size_t N=argc-1, M=N;;
  double input[N];
  for(size_t i=0; i<N; i++)
    input[i] = strtod(argv[i+1],0);

  for(int w=0; w<metric_count; w++){
    puts(metric_list[w].name);

    double result[M][2];
    discrete_fourier_transform(&metric_list[w], N, result, input);
    for(size_t i=0; i<M; i++)
      pnum(result[i][0]);
    puts("");
    for(size_t i=0; i<M; i++)
      pnum(result[i][1]);
    puts("");
    for(size_t i=0; i<N; i++)
      pnum(input[i]);
    puts("");
    double reversed[N];
    inverse_discrete_fourier_transform(&metric_list[w], N, reversed, result);
    for(size_t i=0; i<N; i++)
      pnum(reversed[i]);
    puts("");

    puts("");
  }
}