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

/* 
Demonstration of arrays in C.

Compiles on a UNIX system using 

  gcc -lm -Wall -o C-examples-3 C-examples-3.c

Execute from the shell using "C-examples-3" or "./C-examples-3".
*/


/*****************************************************************************
1-dimensional arrays (vectors)
******************************************************************************/


/* 
Returns the sum of the values in X, where X is an array with length n.
*/
double sum(const double* X,
	   const int n)
{
  int i;
  double S = 0;
  for (i=0; i<n; ++i) S += X[i];
  return S;
}


/* 
Computes the partial sums of X in place, where X is an array with
length n.
*/
void psum(double* X,
	  const int n)
{
  int i;
  for (i=1; i<n; ++i) X[i] += X[i-1];
}


/* 
Computes the first differences of X in place.  Note that the first
component of the first difference vector is undefined.
*/
void diff1(double* X,
	   const int n)
{
  int i;
  for (i=n-1; i>0; --i) X[i] -= X[i-1];
}


/* 
Returns the p-norm length of X, where X has length n.
*/
double pnorm(const double* X,
	     const int n,
	     const double p)
{
  int i;
  double L = 0;
  for (i=0; i<n; ++i) L += exp(p * log(fabs(X[i])));
  return exp(log(L) / p);
}


/* 
Returns the largest value in X, where X has length n. 
*/
double max(const double* X,
	   const int n)
{
  int i;
  double M = X[0];
  for (i=1; i<n; ++i) M = (X[i] > M) ? X[i] : M;
  return M;
}


/* 
Returns the index of the smallest value in X, where X has length n.
*/
int minloc(const double* X,
	   const int n)
{
  int i, L=0;
  for (i=1; i<n; ++i)
    {
      if (X[i] < X[L]) L = i;
    }

  return L;
}


/* 
Returns the inner product of X and Y, where X and Y both have length
n.
*/
double dot(const double* X,
	   const double* Y,
	   const int n)
{
  int i;
  double D = 0;
  for (i=0; i<n; ++i) D += X[i] * Y[i];
  return D;
}


/* 
The "seive of Erasthones" prints all prime numbers between 2 and n.
*/
void seive(const int n)
{
  int i, j;
  
  // Initialize A to hold the values 2, 3, 4, ...
  int* A = (int*)malloc((n-1)*sizeof(int));

  for (i=0; i<n-1; ++i) A[i] = i+2;
  
  for (i=0; i<n-1; ++i)
    {
      /* A[i] has been crossed out, hence so have its multiples. */
      if (A[i] == 0) continue;

      /* Cross out the multiples of A[i]. */
      j = 2;
      while (j*A[i]-2 < n)
	{
	  A[j*A[i]-2] = 0;
	  ++j;
	}
    }
  
  /* Print out the results. */
  for (i=0; i<n-1; ++i) if (A[i] != 0) printf("%4d\n", A[i]);
  
  /* Free up the space that held the sieve. */
  free(A);
}



/*****************************************************************************
2-dimensional arrays (matrices)
******************************************************************************/

/* 
Compute the column sums of A and puts them in S.  A is an m x n array
packed row-wise and S is an n-vector.  Note that we may place the i
loop and the j loop in either order.  The order used here allows the
elements of A to be accessed in order, which in principal should be
slightly faster than doing it the other way.
*/
void column_sum(const double* A,
		const int m,
		const int n,
		double* S)
{
  int i, j;
  memset(S, 0, n * sizeof(double));

  for (i=0; i<m; ++i)
    {
      for (j=0; j<n; ++j) S[j] += A[i*n + j];
    }
}


/* 
Computes the matrix-vector product AX and puts the result in Y.  A is
an m x n array stored row-wise, X is an n-vector, and Y is an
m-vector.
*/
void mvec(const double* A,
	  const double* X,
	  double* Y,
	  const int m,
	  const int n)
{
  int i, j;
  memset(Y, 0, m * sizeof(double));

  for (i=0; i<m; ++i)
    {
      for (j=0; j<n; ++j) Y[i] += A[i*n + j] * X[j];
    }
}


/* 
Computes the matrix-vector product AX and puts the result in Y.  A is
a symmetric matrix with only the lower triangle stored.
*/
void mvecs(const double* A,
	   const double* X,
	   double* Y,
	   const int m)
{
  int  i, j;
  memset(Y, 0, m * sizeof(double));

  for (i=0; i<m; ++i)
    {
      for (j=0; j<i; ++j)
	{
	  Y[i] += A[i*m + j] * X[j];
	  Y[j] += A[i*m + j] * X[i];
	}
      Y[i] += A[i*m + i] * X[i];
    }
}


/* 
Computes the matrix-vector product AX and puts the result in Y.  A is
a symmetric matrix stored in packed form (lower-triangle row-wise or
equivalently upper triangle column-wise).
*/
void mvecp(const double* A,
	   const double* X,
	   double* Y,
	   const int m)
{
  int ii=0, i, j;
  memset(Y, 0, m * sizeof(double));

  for (i=0; i<m; ++i)
    {
      for (j=0; j<i; ++j)
	{
	  Y[i] += A[ii] * X[j];
	  Y[j] += A[ii++] * X[i];
	}
      Y[i] += A[ii++] * X[i];
    }
}


/* 
Computes the matrix-vector product AX and puts the result in Y. A is a
representation of a matrix M in sparse form, so that M(I(k),J(k)) =
A(k), and M(i,j) = 0 if (i,j) != (I(k),J(k)) for some k < q.
*/
void mvecsp(const double* A,
	    const int* I,
	    const int* J,
	    const double* X,
	    double* Y,
	    const int q,
	    const int m,
	    const int n)
{
  int ii;
  memset(Y, 0, m * sizeof(double));
  
  for (ii=0; ii<q; ++ii) Y[I[ii]] += A[ii] * X[J[ii]];
}


int main(int argc, 
	 char** argv)
{
  double *M, *W, *X, *Y, *Z;
  int i, j;

  /* Create a 20-vector and put a quadratic sequence into it. */
  X = (double*)malloc(20*sizeof(double));
  for (i=0; i<20; ++i) X[i] = (i-7) * (i-7) + 3;

  /* Create another 20-vector and put a linear sequence into it. */
  Y = (double*)malloc(20*sizeof(double));
  for (i=0; i<20; ++i) Y[i] = -11 + 4 * i;

  /* Print X and Y in a two-column table. */
  printf("      X               Y\n");
  for (i=0; i<20; ++i) printf("%8.1f\t%8.1f\n", X[i], Y[i]);

  /* Run the non-destructive vector functions. */
  printf("\nsum(X, 20) = %f\n", sum(X, 20));
  printf("dot(X, Y, 20) = %f\n", dot(X, Y, 20));
  printf("pnorm(X, 20, 2) = %f\n", pnorm(X, 20, 2));
  printf("pnorm(X, 20, 7) = %f\n", pnorm(X, 20, 7));
  printf("max(X, 20) = %f\n", max(X, 20));
  printf("minloc(X, 20) = %d\n", minloc(X, 20));

  /* Print prime numbers between 2 and 100. */
  printf("\nThe prime numbers between 2 and 100 are:\n");
  seive(100);

  /* Run the partial sum function.  We want to preserve X so we'll
     make a copy of it to pass into psum(), which is destructive. */
  Z = (double*)malloc(20*sizeof(double));
  memcpy(Z, X, 20 * sizeof(double));
  psum(Z, 20);
  printf("\n        X           psum(X)\n");
  for (i=0; i<20; ++i) printf("%10.1f\t%10.1f\n", X[i], Z[i]);

  /* Run the diff1() function on a copy of X and on a copy of Y. */
  memcpy(Z, X, 20 * sizeof(double));
  W = (double*)malloc(20*sizeof(double));
  memcpy(W, Y, 20 * sizeof(double));
  diff1(Z, 20);
  diff1(W, 20);
  printf("\n        X          diff1(X)             Y          diff1(Y)\n");
  for (i=0; i<20; ++i)
    {
      printf("%10.1f\t%10.1f\t%10.1f\t%10.1f\n", X[i], Z[i], Y[i], W[i]);
    }
  
  /* Construct a 5 x 5 Toeplitz array with quadratically-decreasing
     elements. */
  M = (double*)malloc(25*sizeof(double));
  for (i=0; i<5; ++i)
    {
      for (j=0; j<5; ++j) M[5*i + j] = 1 / (double)(1 + (i-j) * (i-j));
    }

  /* Compute the product MX and place the result into Y. */
  mvec(M, X, Y, 5, 5);
  printf("\n\n");
  for (i=0; i<5; ++i)
    {
      for (j=0; j<5; ++j) printf("%6.2f", M[i*5 + j]);
      printf("\t%s\t", (i == 2) ? "  *" : "   ");
      printf("%5.2f", X[i]);
      printf("\t%s\t", (i == 2) ? "  =" : "   ");
      printf("%5.2f\n", Y[i]);
    }

  return 0;
}