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

#include "hslift.c"

float
l2norm(float x[], int N)
{
  int i;
  double sum;
  sum=0;
  for(i=0; i<N; i++) sum += x[i]*x[i];
  return (float)sqrt(sum);
}

float
l2diff(float x[], float y[], int N)
{
  int i;
  double diff, sum;
  sum=0;
  for(i=0; i<N; i++) {
    diff = x[i]-y[i];
    sum += diff*diff;
  }
  return (float)sqrt(sum);
}


void
make_unit( float u[], int N ) /* Fill v[] with a random unit vector */
{
  int i; float norm;
  for(i=0; i<N; i++)
    u[i] = 1.0 - 2.0*rand()/(float)RAND_MAX; /* in [-1,1] */
  norm = 1.0/l2norm(u,N);
  for(i=0; i<N; i++) u[i] *= norm;
  return;
}

void
copy( float v[], float u[], int N ) /* Copy u[] into v[] */
{
  int i;
  for(i=0; i<N; i++) v[i] = u[i];
  return;
}

int
main(void)
{
  float u[1000]={0}, v[1000]={0};
  int i, N, J;

  N=145;			/* Not a power of 2 */
  J = 8;			/* ceiling(log_2(N)) */

  printf("%d point Haar transform by half-sample symmetric lifting.\n", N);

  make_unit(u,N);  copy(v,u,N);

  hsldht(u, N, 1, J);

  printf("||u||-||Hu|| = %10.2e;  epsilon_f = %10.2e\n",
	 l2norm(v,N)-l2norm(u,N), FLT_EPSILON);

  hslidht(u, N, 1, J);

  printf("||u-H*Hu|| = %10.2e;  epsilon_f = %10.2e\n",
	 l2diff(v,u,N), FLT_EPSILON);

  return 0;
}