fft_test.cpp

/******************************************************************************
* MODULE     : fft_test.cpp
* DESCRIPTION: Test FFT
* COPYRIGHT  : (C) 2005  Joris van der Hoeven
*******************************************************************************
* 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 <typeinfo>
#include <numerix/modular.hpp>
#include <algebramix/fft_naive.hpp>
#include <algebramix/fft_simd.hpp>
#include <algebramix/fft_blocks.hpp>
#include <algebramix/fft_threads.hpp>

using namespace mmx;

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[bit_mirror(i,n)], sum));
  }
}

template<typename C, typename FFTER> void
test_fft () {
  nat n,j;
  C* a, * b, *c;
  FFTER dummy (2);
  mmout << "Testing " << typeid (dummy).name () << "\n";
  nat k= primitive_root_max_order<C> (2);
  nat mx= (nat) 1 << 15;
  mx= k == 0 ? mx : min (mx, primitive_root_max_order<C> (2));
  
  for (n=2; n <= mx;  n<<=1) {
    mmout << "\rwith " << n << " entries" << flush_now;
    typedef implementation<vector_linear,vector_naive> NVec;
    FFTER ffter (n);
    a= make<C> (n);
    b= mmx_new<C> (n); NVec::copy (b, a, n);
    c= mmx_new<C> (n); NVec::copy (c, a, n);

    ffter.direct_transform (b);
    if (n<3000) check (b, a, n);

    for (nat i=0; i < ffter.depth; i++)
      ffter.dfft (c, 1, 0, ffter.depth, i, i+1);
    if (n<3000) check (c, a, n);

    ffter.inverse_transform (b);
    for (j=0; j<n; ++j)
      assert (approx_eq (a[j], b[j]));

    for (nat i=ffter.depth-1; (int) i >= 0 ; i--)
      ffter.ifft (c, 1, 0, ffter.depth, i, i+1);
    for (j=0; j<n; ++j)
      assert (approx_eq (a[j], c[j]/C(n)));

    mmx_delete<C> (a, n);
    mmx_delete<C> (b, n);
    mmx_delete<C> (c, n);
  }
  mmout << "\r";
}

template<typename C> void
test_all_ffts () {

  test_fft<C, fft_naive_transformer<C> >();

  test_fft<C, fft_blocks_transformer<C> >();
  test_fft<C, fft_blocks_transformer<C,
    fft_naive_transformer<C>, 1, 0, 2, 2> >();
  test_fft<C, fft_threads_transformer<C> >();

  test_fft<C, fft_simd_transformer<C> >();
  test_fft<C, fft_simd_transformer<C, fft_blocks_transformer<C> > >();

  test_fft<C, fft_threads_transformer<C,
    fft_blocks_transformer<C, fft_simd_transformer<C>,
      8, 5, 10, 16>, 13 > > ();

}

int
main () {
  try {
    test_all_ffts<complex<double> > ();
    test_all_ffts<modular<modulus<nat>, FFTP> > ();
    test_all_ffts<modular<modulus<unsigned char>,
                          modular_fixed<nat,251> > > ();
    test_all_ffts<modular<modulus<unsigned char>,
                          modular_fixed<nat,193> > > ();
    test_all_ffts<modular<modulus<int>, modular_fixed<int,(1<<16)+1> > > ();
  }
  catch (exception& e) {
    mmout << e << "\n";
    return -1;
  }
  return 0;
}