#include <stdio.h>

class SquareMatrix
{
  int len;

public:
  SquareMatrix (int new_len)
  {
    len = new_len;
  }
  int getLen () const
  {
    return len;
  }
  virtual double getElement (int i, int j) const = 0;
  virtual void setElement (int i, int j, double d) = 0;
};

class DenseMatrix : public SquareMatrix
{
  double *data;

public:
  DenseMatrix (int new_len) : SquareMatrix (new_len)
  {
    data = new double[new_len * new_len];
  }

  double getElement (int i, int j) const;
  void setElement (int i, int j, double d);
  void fill (double d);

  virtual ~DenseMatrix ()
  {
    delete[] data;
  }


};

double DenseMatrix::getElement (int i, int j) const
{
  int len = getLen ();
  return data[i*len+j];
}

void DenseMatrix::setElement (int i, int j, double d)
{
  int len = getLen ();
  data[i*len+j] = d;
}

void DenseMatrix::fill (double d)
{
  int len = getLen ();
  int i, j;

  for (i = 0; i < len; i++)
    {
      for (j = 0; j < len; j++)
	setElement (i, j, d);
      d += 0.5;
    }
}

static void multiply_matrices (SquareMatrix *res, SquareMatrix *a, SquareMatrix *b)
{
  int i, j, len = res->getLen ();

  for (i = 0; i < len; i++)
    for (j = 0; j < len; j++)
      {
	int k;
	double d = 0;

	for (k = 0; k < len; k++)
	  d += a->getElement (i, k) * b->getElement (k, j);

	res->setElement (i, j, d);
      }
}

#define N 1000

int main (int argc, char *argv[])
{
  DenseMatrix r(N), a(N), b(N);

  a.fill(10);
  b.fill(-4);

  multiply_matrices (&r, &a, &b);

  printf ("Some value: %e\n", r.getElement (N/2, N/3));
  return 0;
}