#include <stdio.h>
#include <stdlib.h>

int compute(int N) {

  int i_0, i_1;
  // claim space
  double **A = (double **)malloc(sizeof(double *) * N);
  for (i_0 = 0; i_0 < N; i_0++) {
    A[i_0] = (double *)malloc(sizeof(double) * N);
  }
  double *x = (double *)malloc(sizeof(double) * N);
  double *y = (double *)malloc(sizeof(double) * N);
  double *y_sol = (double *)malloc(sizeof(double) * N);

  // fill matrix and vector
  for (i_0 = 0; i_0 < N; i_0++) {
    x[i_0] = i_0 + 1;
    y_sol[i_0] = (i_0 + 1) * N;
    for (i_1 = 0; i_1 < N; i_1++) {
      A[i_0][i_1] = (i_0 + 1.0) / (i_1 + 1.0);
    }
  }

  // for (i_0 = 0; i_0 < N; i_0++){
  // 	printf("[ ");
  // 	for(i_1 = 0; i_1 < N; i_1++){
  // 		printf("%.3f ",A[i_0][i_1]);
  // 	}
  // 	printf("]");
  // 	printf("[ %.1f ] = [ %.1f ][ %.1f ]\n",x[i_0],y[i_0],y_sol[i_0]);
  // }

  // compute matrix-vector product
  for (i_0 = 0; i_0 < N; i_0++) {
    y[i_0] = 0;
    for (i_1 = 0; i_1 < N; i_1++) {
      y[i_0] += A[i_0][i_1] * x[i_1];
    }
  }

  // compare result to solution
  double sum = 0.0;
  for (i_0 = 0; i_0 < N; i_0++) {
    sum += (y[i_0] - y_sol[i_0]) * (y[i_0] - y_sol[i_0]);
  }

  return sum;
}

int main(int argc, char *argv[]) {

  int N;
  if (argc != 2) {
    N = 1000;
  } else {
    N = atoi(argv[1]); // Convert the argument to an integer
  }

  double err = compute(N);
  printf("Error for N = %d: %.8f\n", N, err);

  return 0;
}