include/algebramix/fft_roots.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_roots.hpp
* DESCRIPTION: primitive roots for diadic and triadic FFTs
* 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.
******************************************************************************/

#ifndef __MMX__FFT_ROOTS__HPP
#define __MMX__FFT_ROOTS__HPP
#include <basix/int.hpp>
#include <basix/vector_naive.hpp>
#include <numerix/modular.hpp>
#include <numerix/modular_int.hpp>
#include <numerix/complex_double.hpp>

namespace mmx {

/******************************************************************************
* Helper functions on bits
******************************************************************************/

// diadic
static nat
bit_mirror (nat i, nat n) {
  if (n == 1) return i;
  else return (bit_mirror (i & ((n>>1) -1), n>>1) << 1) + i / (n>>1);
}

// triadic
static nat
digit_mirror_triadic (nat i, nat n) {
  if (n == 1) return i;
  else return digit_mirror_triadic (i % (n / 3), n / 3) * 3 + i / (n / 3);
}

/******************************************************************************
* Primitive roots of unity
******************************************************************************/

template<typename C>
struct primitive_root_helper {

  // max order for b-adic fft
  static inline nat
  max_order (nat b) { return 0; } // means infinity
    
  static inline C
  op (nat n, nat i, const format<C>& fm) {
    if (n == 1 || i == 0) return promote (1, fm);
    else if (n == 2) return promote (-1, fm);
    else {
      C z= promote (1, fm);
      set_pi (z);
      z= times_i (z);
      z= (promote ((int) (i+i), z) / promote ((int) n, z)) * z;
      return exp (z);
      //Real_type(C) angle= (2 * i * Pi(Real_type(C))) / n;
      //return gaussian (cos (angle), sin (angle));
    }
  }
};

template<typename C> C
primitive_root (nat n, nat i, const format<C>& fm) {
  return primitive_root_helper<C>::op (n, i, fm);
}

template<typename C> nat
primitive_root_max_order (nat b) {
  return primitive_root_helper<C>::max_order (b);
}

// double
#ifdef ACCURATE_ROOTS
template<>
struct primitive_root_helper<complex<double> > {
  static inline nat
  max_order (nat b) { return 0; }

  static complex<double> op (nat n, nat i, const format<complex<double> >&);
};
#endif // ACCURATE_ROOTS

// modular native integers
template<typename I> nat
primitive_root_max_order_int (nat b, I p, I& m) {
  // p must be nonnegative
  nat k= 0;
  m= (nat) (p-1);
  while (m % b == 0) { k++; m /= b; }
  return ((nat) (p-1)) / m;
}

template<typename I> I
primitive_root_max_int (nat b, I p, nat& k, I& m) {
  // root of maximal order k for radix b modulo p
  typedef modulus<I, modulus_int_preinverse<8*sizeof(I)> > Modulus;
  typedef modular<Modulus> Modular;
  k= primitive_root_max_order_int (b, p, m);
  if (k == 1) return I (1);
  Modulus tmp= Modular::get_modulus ();
  Modular::set_modulus (p);
  Modular v;
  for (I x = 1; x < p; x++) {
    v = binpow (Modular (x), (nat) m);
    if (v == 1) continue;
    if (binpow (v, k / b) != 1) break;
  }
  Modular::set_modulus (tmp);
  return * v;
}

template<typename I, typename V, typename W>
struct primitive_root_helper_modular_int {
  typedef modular<modulus<I,V>,W> C;
  
  static inline nat
  max_order (nat b) {
    I m, p= * C::get_modulus ();
    return primitive_root_max_order_int (b, p, m);
  }
    
  static inline C
  op (nat n, nat i, const format<C>& fm) {
    (void) fm; /*FIXME*/ // the format should be used
    nat b, k;
    if (n == next_power_of_two (n)) b = 2;
    else if (n == next_power_of_three (n)) b = 3;
    else ERROR ("only radices 2 and 3 are implemented");
    I m, p= * C::get_modulus ();
    C v= primitive_root_max_int (b, p, k, m);
    VERIFY (n <= k, "maximal order exceeded");
    // v is now a primitive root of maximal order k
    return binpow (v, i * (k / n));
  }
};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<unsigned char,V>,W> > :
    primitive_root_helper_modular_int<unsigned char, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<signed char,V>,W> > :
    primitive_root_helper_modular_int<signed char, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<short unsigned int,V>,W> > :
    primitive_root_helper_modular_int<short unsigned int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<short signed int,V>,W> > :
    primitive_root_helper_modular_int<short signed int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<unsigned int,V>,W> > :
    primitive_root_helper_modular_int<unsigned int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<int,V>,W> > :
    primitive_root_helper_modular_int<int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<long unsigned int,V>,W> > :
    primitive_root_helper_modular_int<long unsigned int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<long int,V>,W> > :
    primitive_root_helper_modular_int<long int, V, W> {};

#ifdef BASIX_HAVE_LONG_LONG_INT
template<typename V, typename W>
struct primitive_root_helper<modular<modulus<long long unsigned int,V>,W> > :
    primitive_root_helper_modular_int<long long unsigned int, V, W> {};

template<typename V, typename W>
struct primitive_root_helper<modular<modulus<long long int,V>,W> > :
    primitive_root_helper_modular_int<long long int, V, W> {};
#endif

// special modulus 3 * 2^30 + 1 
#define FFTP modular_fixed<nat,0x0c0000001> 

template<typename V>
struct primitive_root_helper<modular<modulus<nat,V>,FFTP> > {
  typedef modulus<nat, V> Modulus;
  typedef modular<Modulus,FFTP> Modular;

  static inline nat
  max_order (nat b) {
    if (b == 2) return 1<<30;
    if (b == 3) return 3;
    return 0;
  }

  static inline Modular
  op (nat n, nat i, const format<Modular >&) {
    nat k;
    Modular v;
    if (n == next_power_of_two (n)) {
      k = max_order (2);
      v = 125;
    }
    else if (n == next_power_of_three (n)) {
      k = max_order (3);
      v = 1610563584;
    }
    else ERROR ("only radices 2 and 3 are implemented");
    ASSERT (n <= k, "maximal order exceeded");
    return binpow (v, i * (k / n));
  }
};

// special modulus 7 * 2^26 + 1 
#define FFTP29 modular_fixed<nat,469762049> 

template<typename V>
struct primitive_root_helper<modular<modulus<nat,V>,FFTP29> > {
  typedef modulus<nat, V> Modulus;
  typedef modular<Modulus,FFTP29> Modular;

  static inline nat
  max_order (nat b) {
    if (b == 2) return 1<<26;
    return 0;
  }

  static inline Modular
  op (nat n, nat i, const format<Modular >&) {
    nat k;
    Modular v;
    if (n == next_power_of_two (n)) {
      k = max_order (2);
      v = 2187;
    }
    else { ERROR ("only radix 2 is implemented"); }
    ASSERT (n <= k, "maximal order exceeded");
    return binpow (v, i * (k / n));
  }
};

// special modulus 3 * 2^12 + 1 
#define FFTP14 modular_fixed<unsigned short, 3*(1<<12)+1>

template<typename V>
struct primitive_root_helper<modular<modulus<nat,V>,FFTP14> > {
  typedef modulus<nat, V> Modulus;
  typedef modular<Modulus,FFTP14> Modular;

  static inline nat
  max_order (nat b) {
    if (b == 2) return 1<<12;
    return 0;
  }

  static inline Modular
  op (nat n, nat i, const format<Modular >&) {
    nat k;
    Modular v;
    if (n == next_power_of_two (n)) {
      k = max_order (2);
      v = 1331;
    }
    else { ERROR ("only radix 2 is implemented"); }
    ASSERT (n <= k, "maximal order exceeded");
    return binpow (v, i * (k / n));
  }
};
 
// Bivariate polynomials
template<typename C, typename V> struct polynomial;

template<typename C, typename V>
struct primitive_root_helper<polynomial<C,V> > :
  primitive_root_helper<C> {};

// Polynomials over modulars of polynomials
template<typename C, typename U, typename V, typename W>
struct primitive_root_helper<modular<modulus<polynomial<C,U>,V>,W> > :
  primitive_root_helper<C> {};

/******************************************************************************
* Cache for roots
******************************************************************************/

template<typename V>
struct cached_roots_helper: public V {
  typedef typename V::C C;
  typedef typename V::U U;
  static nat capacity;
  static U*  roots;
  
  static U*
  create_roots (nat n, const format<C>& fm) {
    if (n > capacity) {
      if (capacity != 0) V::destroy_roots (roots, capacity);
      roots= V::create_roots (n, fm);
      capacity= n;
    }
    return roots;
  }

  static void
  destroy_roots (U* u, nat n) {
    (void) u; (void) n;
  }
};

template<typename V> nat cached_roots_helper<V>::capacity= 0;
template<typename V> typename V::U* cached_roots_helper<V>::roots= NULL;

// Root caching is only used for when C can not be modified howsoever
template<typename CC, typename UU=CC, typename SS=CC>
struct roots_helper;

template<typename C>
struct std_roots {
  typedef roots_helper<C> roots_type;
};

template<typename C, typename V>
struct std_roots<vector<C, V> > {
  typedef roots_helper<vector<C, V>, C, C> roots_type;
};

template<typename C, typename V>
struct std_roots<polynomial<C, V> > {
  typedef roots_helper<polynomial<C, V>, C, C> roots_type;
};

STMPL
struct std_roots<complex<double> > {
  typedef complex<double> C;
  typedef cached_roots_helper<roots_helper<C> > roots_type;
};

template<typename M, typename D, D d>
struct std_roots<modular<M, modular_fixed<D, d> > >{
  typedef modular<M, modular_fixed<D, d> > C;
  typedef cached_roots_helper<roots_helper<C> > roots_type;
};

} // namespace mmx
#endif //__MMX__FFT_ROOTS__HPP

---