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

# include "vandermonde_interp_2d.h"

/******************************************************************************/

int triangle_num ( int n )

/******************************************************************************/
/*
  Purpose:

    TRIANGLE_NUM returns the N-th triangular number.

  Discussion:

    The N-th triangular number T(N) is formed by the sum of the first
    N integers:

      T(N) = sum ( 1 <= I <= N ) I

    By convention, T(0) = 0.

    T(N) can be computed quickly by the formula:

      T(N) = ( N * ( N + 1 ) ) / 2

  First Values:

     0
     1
     3
     6
    10
    15
    21
    28
    36
    45
    55

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    07 October 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the index of the desired number, 
    which must be at least 0.

    Output, int TRIANGLE_NUM, the N-th triangular number.
*/
{
  int value;

  value = ( n * ( n + 1 ) ) / 2;

  return value;
}
/******************************************************************************/

double *vandermonde_interp_2d_matrix ( int n, int m, double x[], double y[] )

/******************************************************************************/
/*
  Purpose:

    VANDERMONDE_INTERP_2D_MATRIX computes a Vandermonde 2D interpolation matrix.

  Discussion:

    We assume the approximating function has the form of a polynomial
    in X and Y of total degree M.

      p(x,y) = c00 
             + c10 * x                + c01 * y
             + c20 * x^2   + c11 * xy + c02 * y^2
             + ...
             + cm0 * x^(m) + ...      + c0m * y^m.

    If we let T(K) = the K-th triangular number 
            = sum ( 1 <= I <= K ) I
    then the number of coefficients in the above polynomial is T(M+1).

    We have n data locations (x(i),y(i)) and values z(i) to approximate:

      p(x(i),y(i)) = z(i)

    and we assume that N = T(M+1).

    This can be cast as an NxN linear system for the polynomial
    coefficients:

      [ 1 x1 y1  x1^2 ... y1^m ] [ c00 ] = [  z1 ]
      [ 1 x2 y2  x2^2 ... y2^m ] [ c10 ] = [  z2 ]
      [ 1 x3 y3  x3^2 ... y3^m ] [ c01 ] = [  z3 ]
      [ ...................... ] [ ... ] = [ ... ]
      [ 1 xn yn  xn^2 ... yn^m ] [ c0n ] = [  zn ]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    21 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of data points.  It is necessary 
    that N = T(M+1), where T(K) is the K-th triangular number.

    Input, int M, the degree of the polynomial.

    Input, double X[N], Y[N], the data locations.

    Output, double VANDERMONDE_INTERP_2D_MATRIX[N*N], the Vandermonde matrix for X.
*/
{
  double *a;
  int ex;
  int ey;
  int i;
  int j;
  int s;
  int tmp1;

  tmp1 = triangle_num ( m + 1 );

  if ( n != tmp1 )
  {
    printf ( "\n" );
    printf ( "VANDERMONDE_INTERP_2D_MATRIX - Fatal error!\n" );
    printf ( "  For interpolation, we need N = T(M+1).\n" );
    printf ( "  But we have N = %d\n", n );
    printf ( "  M = %d\n", m );
    printf ( "  and T(M+1) = %d\n", tmp1 );
    exit ( 1 );
  }

  a = ( double * ) malloc ( n * n * sizeof ( double ) );
  j = 0;

  for ( s = 0; s <= m; s++ )
  {
    for ( ex = s; 0 <= ex; ex-- )
    {
      ey = s - ex;
      for ( i = 0; i < n; i++ )
      {
        a[i+j*n] = pow ( x[i], ex ) * pow ( y[i], ey );
      }
      j = j + 1;
    }
  }
 
  return a;
}