/Users/mourrain/Devel/mmx/algebramix/include/algebramix/fkt_transform.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fkt_transform.hpp
* DESCRIPTION: fast Karatsuba transform
* COPYRIGHT  : (C) 2004  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.
******************************************************************************/

#ifndef __MMX__FKT_TRANSFORM__HPP
#define __MMX__FKT_TRANSFORM__HPP
#include <algebramix/polynomial_naive.hpp>

namespace mmx {
#define TMPL template<class C>

template<typename V>
struct fkt_package {
  typedef implementation<polynomial_multiply,V> Pol;
  typedef typename Pol::Vec Vec;
  
TMPL static void
direct_fkt_step (C* dest, nat n) {
  // dest of length n on entry and 3n/2 on exit
  nat h= n >> 1;
  Vec::half_copy (dest + n, dest + 1, h);
  Vec::half_copy (dest + 1, dest + 2, h - 1);
  Pol::add (dest + h, dest, dest + n, h);
}

TMPL static void
direct_fkt (C* dest, nat n2, nat n3) {
  // dest of length n2 = 2^p on entry and n3 = 3^p on exit
  if (n2 == 2) {
    dest[2]= dest[1];
    dest[1] += dest[0];
  }
  else if (n2 == 4) {
    C p= dest[0] + dest[1];
    C q= dest[2] + dest[3];
    dest[6]= dest[2];
    dest[2]= dest[1];
    dest[8]= dest[3];
    dest[1]= p;
    dest[3]= dest[0] + dest[6];
    dest[4]= p+q;
    dest[5]= dest[2] + dest[8];
    dest[7]= q;
  }
  else while (n2 != n3) {
    direct_fkt_step (dest, n2);
    n2= n2 + (n2 >> 1);
  }
}

TMPL static void
inverse_fkt_step (C* dest, nat n, nat l) {
  // dest of length n = 3 k l on entry for some integer k
  // dest of length (2 l + 1) k on exit
  nat t= n/3, k= t/l, ll= (l << 1) + 1;
  Pol::sub (dest + t, dest, t);
  Pol::sub (dest + t, dest + (t << 1), t);
  nat buf_size= aligned_size<C,V> (n);
  C* buf= mmx_new<C> (buf_size);
  Pol::copy (buf, dest, n);
  C* src= buf;
  if (l == 1) {
    Vec::triple_copy (dest, buf, t);
    Vec::triple_copy (dest + 1, buf + t, t);
    Vec::triple_copy (dest + 2, buf + (t << 1), t);
  }
  else while (k != 0) {
    Vec::double_copy (dest, src, l);
    Vec::double_copy (dest + 1, src + t, l);
    Vec::double_add (dest + 2, src + (t << 1), l-1);
    dest[l<<1]= src[(t<<1) + (l-1)];
    dest += ll; src += l; k--;
  }
  mmx_delete<C> (buf, buf_size);
}

TMPL static void
inverse_fkt (C* dest, nat n3, nat n2) {
  // dest of length n3 = 3^p on entry and 2 n2 - 1 on exit with n2 = 2 ^ p
  if (n2 == 2)
    dest[1] -= (dest[0] + dest[2]);
  else if (n2 == 4) {
    dest[3] -= (dest[0] + dest[6]);
    dest[4] -= (dest[1] + dest[7]);
    dest[5] -= (dest[2] + dest[8]);
    C aux= dest[1]; dest[1]= dest[3]; dest[3]= aux;
    aux= dest[2]; dest[2]= dest[6]; dest[6]= aux;
    aux= dest[5]; dest[5]= dest[7]; dest[7]= aux;
    dest[3] -= (dest[0] + dest[6]);
    dest[4] -= (dest[1] + dest[7]);
    dest[5] -= (dest[2] + dest[8]);
    aux= dest[1];
    dest[1]= dest[3];
    dest[3]= dest[4];
    dest[4]= dest[2] + dest[7];
    dest[2]= aux + dest[6];
    dest[6]= dest[8];
  }
  else {
    nat n= n3, l= 1;
    while (n != l) {
      inverse_fkt_step (dest, n, l);
      nat k= n / ((l << 1) + l);
      l= (l << 1) + 1;
      n= k * l;
    }
  }
}

};

#undef TMPL
} // namespace mmx
#endif //__MMX__FKT_TRANSFORM__HPP