algebramix: include/algebramix/fft

00001 
00002 /******************************************************************************
00003 * MODULE     : fft_roots.hpp
00004 * DESCRIPTION: primitive roots for diadic and triadic FFTs
00005 * COPYRIGHT  : (C) 2005  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX__FFT_ROOTS__HPP
00014 #define __MMX__FFT_ROOTS__HPP
00015 #include <basix/int.hpp>
00016 #include <basix/vector_naive.hpp>
00017 #include <numerix/modular.hpp>
00018 #include <numerix/modular_int.hpp>
00019 #include <numerix/complex_double.hpp>
00020 
00021 namespace mmx {
00022 
00023 /******************************************************************************
00024 * Helper functions on bits
00025 ******************************************************************************/
00026 
00027 // diadic
00028 static nat
00029 bit_mirror (nat i, nat n) {
00030   if (n == 1) return i;
00031   else return (bit_mirror (i & ((n>>1) -1), n>>1) << 1) + i / (n>>1);
00032 }
00033 
00034 // triadic
00035 static nat
00036 digit_mirror_triadic (nat i, nat n) {
00037   if (n == 1) return i;
00038   else return digit_mirror_triadic (i % (n / 3), n / 3) * 3 + i / (n / 3);
00039 }
00040 
00041 /******************************************************************************
00042 * Primitive roots of unity
00043 ******************************************************************************/
00044 
00045 template<typename C>
00046 struct primitive_root_helper {
00047 
00048   // max order for b-adic fft
00049   static inline nat
00050   max_order (nat b) { return 0; } // means infinity
00051     
00052   static inline C
00053   op (nat n, nat i, const format<C>& fm) {
00054     if (n == 1 || i == 0) return promote (1, fm);
00055     else if (n == 2) return promote (-1, fm);
00056     else {
00057       C z= promote (1, fm);
00058       set_pi (z);
00059       z= times_i (z);
00060       z= (promote ((int) (i+i), z) / promote ((int) n, z)) * z;
00061       return exp (z);
00062       //Real_type(C) angle= (2 * i * Pi(Real_type(C))) / n;
00063       //return gaussian (cos (angle), sin (angle));
00064     }
00065   }
00066 };
00067 
00068 template<typename C> C
00069 primitive_root (nat n, nat i, const format<C>& fm) {
00070   return primitive_root_helper<C>::op (n, i, fm);
00071 }
00072 
00073 template<typename C> nat
00074 primitive_root_max_order (nat b) {
00075   return primitive_root_helper<C>::max_order (b);
00076 }
00077 
00078 // double
00079 #ifdef ACCURATE_ROOTS
00080 template<>
00081 struct primitive_root_helper<complex<double> > {
00082   static inline nat
00083   max_order (nat b) { return 0; }
00084 
00085   static complex<double> op (nat n, nat i, const format<complex<double> >&);
00086 };
00087 #endif // ACCURATE_ROOTS
00088 
00089 // modular native integers
00090 template<typename I> nat
00091 primitive_root_max_order_int (nat b, I p, I& m) {
00092   // p must be nonnegative
00093   nat k= 0;
00094   m= (nat) (p-1);
00095   while (m % b == 0) { k++; m /= b; }
00096   return ((nat) (p-1)) / m;
00097 }
00098 
00099 template<typename I> I
00100 primitive_root_max_int (nat b, I p, nat& k, I& m) {
00101   // root of maximal order k for radix b modulo p
00102   typedef modulus<I, modulus_int_preinverse<8*sizeof(I)> > Modulus;
00103   typedef modular<Modulus> Modular;
00104   k= primitive_root_max_order_int (b, p, m);
00105   if (k == 1) return I (1);
00106   Modulus tmp= Modular::get_modulus ();
00107   Modular::set_modulus (p);
00108   Modular v;
00109   for (I x = 1; x < p; x++) {
00110     v = binpow (Modular (x), (nat) m);
00111     if (v == 1) continue;
00112     if (binpow (v, k / b) != 1) break;
00113   }
00114   Modular::set_modulus (tmp);
00115   return * v;
00116 }
00117 
00118 template<typename I, typename V, typename W>
00119 struct primitive_root_helper_modular_int {
00120   typedef modular<modulus<I,V>,W> C;
00121   
00122   static inline nat
00123   max_order (nat b) {
00124     I m, p= * C::get_modulus ();
00125     return primitive_root_max_order_int (b, p, m);
00126   }
00127     
00128   static inline C
00129   op (nat n, nat i, const format<C>& fm) {
00130     (void) fm; /*FIXME*/ // the format should be used
00131     nat b, k;
00132     if (n == next_power_of_two (n)) b = 2;
00133     else if (n == next_power_of_three (n)) b = 3;
00134     else ERROR ("only radices 2 and 3 are implemented");
00135     I m, p= * C::get_modulus ();
00136     C v= primitive_root_max_int (b, p, k, m);
00137     VERIFY (n <= k, "maximal order exceeded");
00138     // v is now a primitive root of maximal order k
00139     return binpow (v, i * (k / n));
00140   }
00141 };
00142 
00143 template<typename V, typename W>
00144 struct primitive_root_helper<modular<modulus<unsigned char,V>,W> > :
00145     primitive_root_helper_modular_int<unsigned char, V, W> {};
00146 
00147 template<typename V, typename W>
00148 struct primitive_root_helper<modular<modulus<signed char,V>,W> > :
00149     primitive_root_helper_modular_int<signed char, V, W> {};
00150 
00151 template<typename V, typename W>
00152 struct primitive_root_helper<modular<modulus<short unsigned int,V>,W> > :
00153     primitive_root_helper_modular_int<short unsigned int, V, W> {};
00154 
00155 template<typename V, typename W>
00156 struct primitive_root_helper<modular<modulus<short signed int,V>,W> > :
00157     primitive_root_helper_modular_int<short signed int, V, W> {};
00158 
00159 template<typename V, typename W>
00160 struct primitive_root_helper<modular<modulus<unsigned int,V>,W> > :
00161     primitive_root_helper_modular_int<unsigned int, V, W> {};
00162 
00163 template<typename V, typename W>
00164 struct primitive_root_helper<modular<modulus<int,V>,W> > :
00165     primitive_root_helper_modular_int<int, V, W> {};
00166 
00167 template<typename V, typename W>
00168 struct primitive_root_helper<modular<modulus<long unsigned int,V>,W> > :
00169     primitive_root_helper_modular_int<long unsigned int, V, W> {};
00170 
00171 template<typename V, typename W>
00172 struct primitive_root_helper<modular<modulus<long int,V>,W> > :
00173     primitive_root_helper_modular_int<long int, V, W> {};
00174 
00175 #ifdef BASIX_HAVE_LONG_LONG_INT
00176 template<typename V, typename W>
00177 struct primitive_root_helper<modular<modulus<long long unsigned int,V>,W> > :
00178     primitive_root_helper_modular_int<long long unsigned int, V, W> {};
00179 
00180 template<typename V, typename W>
00181 struct primitive_root_helper<modular<modulus<long long int,V>,W> > :
00182     primitive_root_helper_modular_int<long long int, V, W> {};
00183 #endif
00184 
00185 // special modulus 3 * 2^30 + 1 
00186 #define FFTP modular_fixed<nat,0x0c0000001> 
00187 
00188 template<typename V>
00189 struct primitive_root_helper<modular<modulus<nat,V>,FFTP> > {
00190   typedef modulus<nat, V> Modulus;
00191   typedef modular<Modulus,FFTP> Modular;
00192 
00193   static inline nat
00194   max_order (nat b) {
00195     if (b == 2) return 1<<30;
00196     if (b == 3) return 3;
00197     return 0;
00198   }
00199 
00200   static inline Modular
00201   op (nat n, nat i, const format<Modular >&) {
00202     nat k;
00203     Modular v;
00204     if (n == next_power_of_two (n)) {
00205       k = max_order (2);
00206       v = 125;
00207     }
00208     else if (n == next_power_of_three (n)) {
00209       k = max_order (3);
00210       v = 1610563584;
00211     }
00212     else ERROR ("only radices 2 and 3 are implemented");
00213     ASSERT (n <= k, "maximal order exceeded");
00214     return binpow (v, i * (k / n));
00215   }
00216 };
00217 
00218 // Bivariate polynomials
00219 template<typename C, typename V> struct polynomial;
00220 
00221 template<typename C, typename V>
00222 struct primitive_root_helper<polynomial<C,V> > :
00223   primitive_root_helper<C> {};
00224 
00225 // Polynomials over modulars of polynomials
00226 template<typename C, typename U, typename V, typename W>
00227 struct primitive_root_helper<modular<modulus<polynomial<C,U>,V>,W> > :
00228   primitive_root_helper<C> {};
00229 
00230 /******************************************************************************
00231 * Cache for roots
00232 ******************************************************************************/
00233 
00234 template<typename V>
00235 struct cached_roots_helper: public V {
00236   typedef typename V::C C;
00237   typedef typename V::U U;
00238   static nat capacity;
00239   static U*  roots;
00240   
00241   static U*
00242   create_roots (nat n, const format<C>& fm) {
00243     if (n > capacity) {
00244       if (capacity != 0) V::destroy_roots (roots, capacity);
00245       roots= V::create_roots (n, fm);
00246       capacity= n;
00247     }
00248     return roots;
00249   }
00250 
00251   static void
00252   destroy_roots (U* u, nat n) {
00253     (void) u; (void) n;
00254   }
00255 };
00256 
00257 template<typename V> nat cached_roots_helper<V>::capacity= 0;
00258 template<typename V> typename V::U* cached_roots_helper<V>::roots= NULL;
00259 
00260 // Root caching is only used for when C can not be modified howsoever
00261 template<typename CC, typename UU=CC, typename SS=CC>
00262 struct roots_helper;
00263 
00264 template<typename C>
00265 struct std_roots {
00266   typedef roots_helper<C> roots_type;
00267 };
00268 
00269 template<typename C, typename V>
00270 struct std_roots<polynomial<C, V> > {
00271   typedef roots_helper<polynomial<C, V>, C, C> roots_type;
00272 };
00273 
00274 STMPL
00275 struct std_roots<complex<double> > {
00276   typedef complex<double> C;
00277   typedef cached_roots_helper<roots_helper<C> > roots_type;
00278 };
00279 
00280 template<typename M, typename D, D d>
00281 struct std_roots<modular<M, modular_fixed<D, d> > >{
00282   typedef modular<M, modular_fixed<D, d> > C;
00283   typedef cached_roots_helper<roots_helper<C> > roots_type;
00284 };
00285 
00286 } // namespace mmx
00287 #endif //__MMX__FFT_ROOTS__HPP
include/algebramix/fft_roots.hpp
