fft_triadic_test.cpp

/******************************************************************************
* MODULE     : fft_triadic_test.cpp
* DESCRIPTION: Test triadic FFT
* COPYRIGHT  : (C) 2008  Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <basix/vector_naive.hpp>
#include <basix/double.hpp>
#include <numerix/complex.hpp>
#include <numerix/modular_int.hpp>
#include <algebramix/fft_triadic_naive.hpp>
#include <algebramix/fft_triadic_threads.hpp>

using namespace mmx;

#define K complex<double>

template<typename C> bool
approx_eq (const C& x, const C& y) {
  return x == y; }

template<> bool
approx_eq (const complex<double>& x, const complex<double>& y) {
  return abs (x-y) < max (1.0e-2 * max (abs (x), abs (y)), 1.0e-2); }

template<typename C> C*
make (nat n) {
  C* a= mmx_new<C> (n);
  for (nat i=0; i<n; i++)
    a[i]= i*i+1;
  return a;
}

template<typename C> void
check (C* a, C* b, nat n) {
  C u= primitive_root<C> (n, 1);
  C v= 1;
  for (nat i=0; i<n; i++, v *= u) {
    C sum= 0;
    for (int j=n-1; j>=0; j--)
      sum= b[j] + v * sum;
    assert (approx_eq (a[digit_mirror_triadic(i,n)], sum));
  }
}

void
test_fft_triadic () {
  typedef implementation<vector_linear,vector_naive> NVec;
  int n,j;
  K* a;
  K* b;
  for (n=3; n < binpow(3, 15); n *=3) {
    mmout << "Testing triadic FFT on " << n << " entries\r" << flush_now;
    //fft_triadic_naive_transformer<K> ffter (n);
    fft_triadic_threads_transformer<K> ffter (n);
    a= make<K> (n);
    b= mmx_new<K> (n); NVec::copy (b, a, n);
    //NVec::print (mmout, a, n); mmout << "\n";
    ffter.direct_transform_triadic (a);
    //NVec::print (mmout, a, n); mmout << "\n";
    if (n<3000) check (a, b, n);
    ffter.inverse_transform_triadic (a);
    //NVec::print (mmout, a, n); mmout << "\n";
    for(j=0; j<n; ++j)
      assert (approx_eq (a[j], b[j]));
    mmx_delete<K> (a, n);
    mmx_delete<K> (b, n);
  }
  mmout << "\n";
}

int
main () {
  try {
    test_fft_triadic();
  }
  catch (exception& e) {
    mmout << e << "\n";
    return -1;
  }
  return 0;
}