#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fftw3.h>

static const int nx = 8;
static const int ny = 8;

int xy_col_major(int x, int y, int number_of_rows) { return y + x * number_of_rows; }

// Compute the mathematically correct result to pow(-1.0, n) for integer n.  The
// function name "m1_pow_n" abbreviates "minus 1 to the power of n".
double m1_pow_n(int n) { return (n % 2)? -1.0: +1.0; }
double squared(double x) { return x * x; }

struct chebyshev_matrix { double* D; double* x; };
chebyshev_matrix cheb2(int const n)
{
  int const dim = n + 1; // dimension (width & height) of the resulting matrix.

  double* D  = fftw_alloc_real(squared(dim));
  double* xv = fftw_alloc_real(dim); // Renamed x to xv to avoid name conflicts with row index x.
  double* c  = fftw_alloc_real(dim);
  double* X  = fftw_alloc_real(squared(dim));
  double* dX = fftw_alloc_real(squared(dim));
  if (!D or !xv or !c or !X or !dX) std::abort();
  
  if (n == 0) { D[0] = 0.0; xv[0] = 1.0; return { D, xv }; }
  
  for (int i = 0; i < dim; ++i) xv[i] = std::cos(i * M_PI / n);
  for (int i = 0; i < dim; ++i) c[i] = m1_pow_n(i);
  c[0] = 2.0 * m1_pow_n(0);
  c[n] = 2.0 * m1_pow_n(n);

  for (int x = 0; x < dim; ++x)
    for (int y = 0; y < dim; ++y)
      X[xy_col_major(x, y, dim)] = xv[y];

  for (int x = 0; x < dim; ++x)
    for (int y = 0; y < dim; ++y)
      dX[xy_col_major(x, y, dim)] = X[xy_col_major(x, y, dim)] - X[xy_col_major(y, x, dim)];

  for (int x = 0; x < dim; ++x)
    for (int y = 0; y < dim; ++y)
    {
      D[xy_col_major(x, y, dim)] = (c[y] * (1.0 / c[x]));
      D[xy_col_major(x, y, dim)] /= (dX[xy_col_major(x, y, dim)] + (x == y));
    }

  for (int y = 0; y < dim; ++y)
  {
    double sum_of_row = 0.0;
    for (int x = 0; x < dim; ++x)
      sum_of_row += D[xy_col_major(x, y, dim)];
    D[xy_col_major(y, y, dim)] -= sum_of_row;
  }

  fftw_free(dX);
  fftw_free(X);
  fftw_free(c);
  
  return { D, xv };
}

void print(double* m, int width, int height)
{
  std::cout << std::fixed << std::setprecision(4);
  for (int y = 0; y < height; ++y)
  {
    for (int x = 0; x < width; ++x)
      std::cout << std::setw(9) << m[xy_col_major(x, y, width)] << ' ';
    std::cout << '\n';
  }

  std::cout << '\n';
}

int main()
{	
  for (int n = 0; n < 6; ++n)
  {
    chebyshev_matrix const result = cheb2(n);
    double* D = result.D;
    double* x = result.x;
    print(D, n + 1, n + 1);
    fftw_free(x);
    fftw_free(D);
  }
}

#include <complex>
#include <iostream>
#include <iomanip>

#include <fftw3.h>

int xy(int x, int y, int number_of_cols) { return x + y * number_of_cols; }

void matrix_product(
  fftw_complex* out,
  fftw_complex const* a, int a_width, int a_height,
  double const*       b, int b_width, int b_height)
{
  for (int x = 0; x < b_width; ++x)
    for (int y = 0; y < a_height; ++y)
    {
      // compute inner product of xth column of b with yth row of a
      
      fftw_complex inner_product {0.0, 0.0};
      for (int i = 0; i < a_width; ++i)
      {
        double const p_real = a[xy(i, y, a_width)][0];
        double const p_imag = a[xy(i, y, a_width)][1];
        double const q = b[xy(x, i, b_width)];
        
        inner_product[0] += p_real * q;
        inner_product[1] += p_imag * q;
      }
      
      out[xy(x, y, b_width)][0] = inner_product[0];
      out[xy(x, y, b_width)][1] = inner_product[1];
    }
}

void print(fftw_complex* m, int width, int height)
{
  std::cout << std::fixed << std::setprecision(2);
  for (int y = 0; y < height; ++y)
  {
    for (int x = 0; x < width; ++x)
      std::cout << std::setw(6) << std::complex{m[xy(x, y, width)][0], m[xy(x, y, width)][1]} <<  ' ';
    std::cout << '\n';
  }

  std::cout << '\n';
}

int main()
{
  fftw_complex a[] {
    {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0},
    {4.0, 0.0}, {5.0, 0.0}, {6.0, 0.0}
  };
  int const a_width = 3;
  int const a_height = 2;

  double b[] {
    7.0,  8.0,
    9.0, 10.0,
    11.0, 12.0
  };
  int const b_width = 2;
  int const b_height = 3;

  int const p_width = b_width;
  int const p_height = a_height;
  fftw_complex p[p_width * p_height];
  matrix_product(p, a, a_width, a_height, b, b_width, b_height);
  print(p, p_width, p_height);
}