/* -*- c++ -*- */
/*
 * Copyright 2002 Free Software Foundation, Inc.
 * 
 * This file is part of GNU Radio
 * 
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3, or (at your option)
 * any later version.
 * 
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cppunit/TestAssert.h>
#include <qa_ccomplex_dotprod_x86.h>
#include <ccomplex_dotprod_x86.h>
#include <string.h>
#include <iostream>
#include <malloc16.h>
#include <sse_debug.h>
#include <cmath>
#include <gr_cpu.h>
#include <random.h>

using std::cerr;

/// Macro for primitive value comparisons
#define assertcomplexEqual(expected0,expected1,actual,delta)	\
  CPPUNIT_ASSERT_DOUBLES_EQUAL (expected0, actual[0], delta);	\
  CPPUNIT_ASSERT_DOUBLES_EQUAL (expected1, actual[1], delta);	


#define	MAX_BLKS	10
#define	FLOATS_PER_BLK	 4

#define	ERR_DELTA	(1e-6)

static float
uniform ()
{
  return 2.0 * ((float) random () / RANDOM_MAX - 0.5);	// uniformly (-1, 1)
}

static void
random_floats (float *buf, unsigned n)
{
  for (unsigned i = 0; i < n; i++)
    buf[i] = rint (uniform () * 32767);
}

static void
zero_floats (float *buf, unsigned n)
{
  for (unsigned i = 0; i < n; i++)
    buf[i] = 0.0;
}

void
ref_ccomplex_dotprod (const float *input,
		      const float *taps, unsigned n_2_ccomplex_blocks,
		      float *result)
{
  float sum0[2] = {0,0};
  float sum1[2] = {0,0};

  do {

    sum0[0] += input[0] * taps[0] - input[1] * taps[1];
    sum0[1] += input[0] * taps[1] + input[1] * taps[0];
    sum1[0] += input[2] * taps[2] - input[3] * taps[3];
    sum1[1] += input[2] * taps[3] + input[3] * taps[2];

    input += 4;
    taps += 4;

  } while (--n_2_ccomplex_blocks != 0);


  result[0] = sum0[0] + sum1[0];
  result[1] = sum0[1] + sum1[1];
}

void 
qa_ccomplex_dotprod_x86::setUp ()
{
  taps = (float *) calloc16Align (MAX_BLKS, 
				  sizeof (float) * FLOATS_PER_BLK);

  input = (float *) calloc16Align (MAX_BLKS,
				   sizeof (float) * FLOATS_PER_BLK);

  if (taps == 0 || input == 0)
    abort ();
}

void
qa_ccomplex_dotprod_x86::tearDown ()
{
  free16Align (taps);
  free16Align (input);
  taps = 0;
  input = 0;
}


void 
qa_ccomplex_dotprod_x86::zb ()	// "zero both"
{
  zero_floats (taps, MAX_BLKS * FLOATS_PER_BLK);
  zero_floats (input, MAX_BLKS * FLOATS_PER_BLK);
}

// 
// t1 
//

void
qa_ccomplex_dotprod_x86::t1_base (ccomplex_dotprod_t ccomplex_dotprod)
{
  float result[2];

  // cerr << "Testing dump_xmm_regs\n";
  // dump_xmm_regs ();

  // test basic cases, 1 block

  zb ();
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (0.0, 0.0, result, ERR_DELTA);

  // vary each input

  zb ();
  input[0] = 1.0;	taps[0] = 1.0; taps[1] = -1.0;
  ccomplex_dotprod (input, taps, 1, result);
  //cerr << result[0] << " " << result[1] << "\n";
  assertcomplexEqual (1.0, -1.0, result, ERR_DELTA);
  
  zb ();
  input[1] = 2.0;	taps[0] = 1.0; taps[1] = -1.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (2.0, 2.0, result, ERR_DELTA);

  zb ();
  input[2] = 3.0;	taps[2] = 1.0; taps[3] = -1.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (3.0, -3.0, result, ERR_DELTA);
  
  zb ();
  input[3] = 4.0;	taps[2] = 1.0; taps[3] = -1.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (4.0, 4.0, result, ERR_DELTA);

  // vary each tap

  zb ();
  input[0] = 1.0;	taps[0] = 0.5; taps[1] = -0.5;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (0.5, -0.5, result, ERR_DELTA);
  
  zb ();
  input[0] = 1.0;	taps[0] = 2.0; taps[1] = -2.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (2.0, -2.0, result, ERR_DELTA);
  
  zb ();
  input[0] = 1.0;	taps[0] = 3.0; taps[1] = -3.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (3.0, -3.0, result, ERR_DELTA);
  
  zb ();
  input[0] = 1.0;	taps[0] = 4.0; taps[1] = -4.0;
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (4.0, -4.0, result, ERR_DELTA);
}

// 
// t2 
//
void
qa_ccomplex_dotprod_x86::t2_base (ccomplex_dotprod_t ccomplex_dotprod)
{
  float result[2];

  zb ();
  input[0] =  1.0; input[1] =  3.0;	taps[0] =  5.0;	taps[1] =  -2.0;

  //1*5-3*-2 =11, 1*-2+3*5=13
 
  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (11.0, 13.0, result, ERR_DELTA);

  //7*5-13*-5 =100, 7*-5+13*5=30

  input[2] =  7.0; input[3] = 13.0;	taps[2] =  5.0;	taps[3] =  -5.0;

  ccomplex_dotprod (input, taps, 1, result);
  assertcomplexEqual (111.0, 43.0, result, ERR_DELTA);

  input[4] = 19; input[5] = -19;	taps[4] = 23.0;	taps[5] = -23.0;

  //19*23--19*-23 =0, 19*-23+-19*23=-874

  ccomplex_dotprod (input, taps, 2, result);
  assertcomplexEqual (111.0, -831.0, result, ERR_DELTA);
  
}

//
// t3
//
void
qa_ccomplex_dotprod_x86::t3_base (ccomplex_dotprod_t ccomplex_dotprod)
{
  srandom (0);	// we want reproducibility

  for (unsigned int i = 0; i < 10; i++){
    random_floats (input, MAX_BLKS * FLOATS_PER_BLK);
    random_floats (taps, MAX_BLKS * FLOATS_PER_BLK);

    // we use a sloppy error margin because on the x86 architecture,
    // our reference implementation is using 80 bit floating point
    // arithmetic, while the SSE version is using 32 bit float point
    // arithmetic.

    float ref[2];
    ref_ccomplex_dotprod (input, taps, MAX_BLKS, ref);
    float calc[2];
    ccomplex_dotprod (input, taps, MAX_BLKS, calc);
    CPPUNIT_ASSERT_DOUBLES_EQUAL (ref[0],
		    	calc[0],
			fabs (ref[0]) * 1e-4);
    CPPUNIT_ASSERT_DOUBLES_EQUAL (ref[1],
		    	calc[1],
			fabs (ref[1]) * 1e-4);
  }
}

void
qa_ccomplex_dotprod_x86::t1_3dnowext ()
{
  if (!gr_cpu::has_3dnowext ()){
    cerr << "No 3DNow!Ext support; not tested\n";
  }
  else
    t1_base (ccomplex_dotprod_3dnowext);
}

void 
qa_ccomplex_dotprod_x86::t2_3dnowext ()
{
  if (!gr_cpu::has_3dnowext ()){
    cerr << "No 3DNow!Ext support; not tested\n";
  }
  else
    t2_base (ccomplex_dotprod_3dnowext);
}

void 
qa_ccomplex_dotprod_x86::t3_3dnowext ()
{
  if (!gr_cpu::has_3dnowext ()){
    cerr << "No 3DNow!Ext support; not tested\n";
  }
  else
    t3_base (ccomplex_dotprod_3dnowext);
}

void
qa_ccomplex_dotprod_x86::t1_3dnow ()
{
  if (!gr_cpu::has_3dnow ()){
    cerr << "No 3DNow! support; not tested\n";
  }
  else
    t1_base (ccomplex_dotprod_3dnow);
}

void 
qa_ccomplex_dotprod_x86::t2_3dnow ()
{
  if (!gr_cpu::has_3dnow ()){
    cerr << "No 3DNow! support; not tested\n";
  }
  else
    t2_base (ccomplex_dotprod_3dnow);
}

void 
qa_ccomplex_dotprod_x86::t3_3dnow ()
{
  if (!gr_cpu::has_3dnow ()){
    cerr << "No 3DNow! support; not tested\n";
  }
  else
    t3_base (ccomplex_dotprod_3dnow);
}

void 
qa_ccomplex_dotprod_x86::t1_sse ()
{
  if (!gr_cpu::has_sse ()){
    cerr << "No SSE support; not tested\n";
  }
  else
    t1_base (ccomplex_dotprod_sse);
}

void 
qa_ccomplex_dotprod_x86::t2_sse ()
{
  if (!gr_cpu::has_sse ()){
    cerr << "No SSE support; not tested\n";
  }
  else
    t2_base (ccomplex_dotprod_sse);
}

void 
qa_ccomplex_dotprod_x86::t3_sse ()
{
  if (!gr_cpu::has_sse ()){
    cerr << "No SSE support; not tested\n";
  }
  else
    t3_base (ccomplex_dotprod_sse);
}