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

# include "ppc_array.h"
# include "ppc_haar.h"
# include "ppc_xmalloc.h"

void haar_transform_vector_forward ( double *v, int n );
void haar_transform_vector_reverse ( double *v, int n );
void haar_transform_matrix_forward ( double **a, int m, int n );
void haar_transform_matrix_reverse ( double **a, int m, int n );

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

void haar_transform_vector ( double *v, int n, int dir )

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

    haar_transform_vector() applies a Haar transform to a vector.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *v: the vector to be transformed in place.

    int n: the number of entries.
    n must be a power of 2;

    int dir: WT_FWD for forward transform, WT_REV for backward transform.

  Output:

    double *v: the transformed vector.
*/
{
  if ( dir == WT_FWD )
  {
    haar_transform_vector_forward ( v, n );
  }
  else if ( dir == WT_REV )
  {
    haar_transform_vector_reverse ( v, n );
  }
  else
  {
    fprintf ( stderr, "haar_transform_vector(): Fatal error!\n" );
    fprintf ( stderr, "  Third input argument illegal.\n" );
    exit ( EXIT_FAILURE );
  }

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

void haar_transform_vector_forward ( double *v, int n )

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

    haar_transform_vector_forward() applies a forward Haar transform to a vector.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *v: the vector to be transformed in place.

    int n: the number of entries.
    n must be a power of 2.

  Output:

    double *v: the transformed vector.
*/
{
  double c = 1.0 / sqrt ( 2.0 );
  int d;
  double h = sqrt ( n );
  int i;
  double x;
  double y;

  for ( i = 0; i < n; i++ )
  {
    v[i] = v[i] / h;
  }

  for ( d = 1; d < n; d = d * 2 )
  {
    for ( i = 0; i < n; i = i + 2 * d )
    {
      x = c * ( v[i] + v[i+d] );
      y = c * ( v[i] - v[i+d] );
      v[i] = x;
      v[i+d] = y;
    }
  }

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

void haar_transform_vector_reverse ( double *v, int n )

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

    haar_transform_vector_reverse() applies a reverse Haar transform to a vector.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *v: the vector to be transformed in place.

    int n: the number of entries.
    n must be a power of 2.

  Output:

    double *v: the transformed vector.
*/
{
  double c = 1.0 / sqrt ( 2.0 );
  int d;
  double h = sqrt ( n );
  int i;
  double x;
  double y;

  for ( d = n / 2; 1 <= d; d = d / 2 )
  {
    for ( i = 0; i < n; i = i + 2 * d )
    {
      x = c * ( v[i] + v[i+d] );
      y = c * ( v[i] - v[i+d] );
      v[i] = x;
      v[i+d] = y;
    }
  }

  for ( i = 0; i < n; i++ )
  {
    v[i] = v[i] * h;
  }

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

void haar_transform_matrix ( double **a, int m, int n, int dir )

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

    haar_transform_matrix() applies a Haar transform to a matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    12 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double **a: the matrix to be transformed in place.

    int m, n: the number of rows and columns.
    m and n must be powers of 2.

    int dir: WT_FWD for forward transform, WT_REV for backward transform.

  Output:

    double **a: the transformed matrix.
*/
{
  if ( dir == WT_FWD )
  {
    haar_transform_matrix_forward ( a, m, n );
  }
  else if ( dir == WT_REV )
  {
    haar_transform_matrix_reverse ( a, m, n );
  }
  else
  {
    fprintf ( stderr, "haar_transform_matrix(): Fatal error!\n" );
    fprintf ( stderr, "  Fourth input argument illegal.\n" );
    exit ( EXIT_FAILURE );
  }

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

void haar_transform_matrix_forward ( double **a, int m, int n )

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

    haar_transform_matrix_forward() applies a forward Haar transform to a matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    12 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *a: the matrix to be transformed in place.

    int m, n: the number of rows and columns.
    m and n must be powers of 2.

  Output:

    double *a: the transformed matrix.
*/
{
  double c = 1.0 / sqrt ( 2.0 );
  int d;
  double h;
  int i;
  int j;
  double x;
  double y;

  for ( i = 0; i < m; i++ )
  {
    haar_transform_vector ( a[i], n, WT_FWD );
  }

  h = sqrt ( m );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      a[i][j] = a[i][j] / h;
    }

    for ( d = 1; d < m; d = d * 2 )
    {
      for ( i = 0; i < m; i = i + 2 * d )
      {
        x = c * ( a[i][j] + a[i+d][j] );
        y = c * ( a[i][j] - a[i+d][j] );
        a[i][j] = x;
        a[i+d][j] = y;
      }
    }
  }

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

void haar_transform_matrix_reverse ( double **a, int m, int n )

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

    haar_transform_matrix_reverse() applies a backward Haar transform to a matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    12 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *a: the matrix to be transformed in place.

    int m, n: the number of rows and columns.
    m and n must be powers of 2.

  Output:

    double *a: the transformed matrix.
*/
{
  double c = 1.0 / sqrt ( 2.0 );
  int d;
  double h;
  int i;
  int j;
  double x;
  double y;

  h = sqrt ( m );

  for ( j = 0; j < n; j++ )
  {
    for ( d = m / 2; 1 <= d; d = d / 2 )
    {
      for ( i = 0; i < m; i = i + 2 * d )
      {
        x = c * ( a[i][j] + a[i+d][j] );
        y = c * ( a[i][j] - a[i+d][j] );
        a[i][j] = x;
        a[i+d][j] = y;
      }
    }
    for ( i = 0; i < m; i++ )
    {
      a[i][j] = a[i][j] * h;
    }
  }

  for ( i = 0; i < m; i++ )
  {
    haar_transform_vector ( a[i], n, WT_REV );
  }

  return;
}