00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 #ifndef __MMX__BASE_DICHO__HPP
00014 #define __MMX__BASE_DICHO__HPP
00015 #include <algebramix/base_naive.hpp>
00016 
00017 namespace mmx {
00018 
00019 
00020 
00021 
00022 
00023 template<typename C>
00024 struct std_base_dicho {
00025   typedef C                           base;
00026   typedef C                           modulus_base;
00027   typedef typename Modulus_variant(C) modulus_base_variant;
00028 };
00029 
00030 
00031 
00032 
00033 
00034 #define Base_dicho_variant(C) base_dicho_variant_helper<C>::BV
00035 
00036 template<typename V>
00037 struct base_dicho : public V {};
00038 
00039 template<typename F, typename V, typename W>
00040 struct implementation<F,V,base_dicho<W> >:
00041   public implementation<F,V,W> {};
00042 
00043 template<typename C>
00044 struct base_dicho_variant_helper {
00045   typedef base_dicho<typename base_naive_variant_helper<C>::BV> BV;
00046 };
00047 
00048 
00049 
00050 
00051 
00052 template<typename V, typename W>
00053 struct implementation<base_transform,V,base_dicho<W> > {
00054 
00055   template<typename C, typename M, typename I> static inline nat
00056   direct (I* c, nat n, const C& a, const M* pow_p) {
00057     if (n == 0 || a == 0) return 0;
00058     if (n == 1) {
00059       c[0]= rem (a, * pow_p[0]);
00060       return 1; }
00061     C q, r;
00062     nat h= (n+1) >> 1, l= log_2 (h); nat m= ((nat) 1) << l;
00063     r= rem (a, * pow_p[l], q);
00064     nat o= direct (c, m, r, pow_p);
00065     if (q == 0) return o;
00066     for (; o < m; o++) c[o]= 0;
00067     return m + direct (c + m, n - m, q, pow_p); }
00068 
00069   template<typename C, typename M, typename I> static inline void
00070   inverse (C& a, const I* c, nat n, const M* pow_p) {
00071     if (n == 0) return;
00072     if (n == 1) { a= c[0]; return; }
00073     if (n == 2) { a= c[0] + C(* pow_p[0]) * c[1]; return; }
00074     nat h= (n+1) >> 1; nat l= log_2 (h); nat m= ((nat) 1) << l; C b;
00075     inverse (a, c    , m    , pow_p);
00076     inverse (b, c + m, n - m, pow_p);
00077     a += * pow_p[l] * b; }
00078 };
00079 
00080 
00081 
00082 
00083 
00084 template<typename V,typename W>
00085 struct implementation<base_transform,V,base_signed<base_dicho<W> > > :
00086     implementation<base_transform,V,base_dicho<W> > {
00087 private:
00088   template<typename C, typename M, typename I> static inline nat
00089   direct (I* c, nat n, const C& a, const M* pow_p, I& carry) {
00090     if (n == 0 || a == 0) { carry= 0; return 0; }
00091     C q, r;
00092     if (n == 1) {
00093       if (a < 0) {
00094         r= rem (- a, * pow_p[0], q);
00095         if (r > * pow_p[0] >> 1) { r -= * pow_p[0]; q += 1; }
00096         neg (r); neg (q);
00097       }
00098       else {
00099         r= rem (a, * pow_p[0], q);
00100         if (r > * pow_p[0] >> 1) { r -= * pow_p[0]; q += 1; }
00101       }
00102       c[0]= as<I> (r); carry= q;
00103       return c[0] == 0 ? 0 : 1;
00104     }
00105     nat h= (n+1) >> 1, l= log_2 (h); nat m= ((nat) 1) << l;
00106     r= rem (abs (a), * pow_p[l], q);
00107     if (a < 0) { neg (r); neg (q); }
00108     nat o= direct (c, m, r, pow_p, carry); q += carry;
00109     if (q == 0) return o;
00110     for (; o < m; o++) c[o]= 0;
00111     return m + direct (c + m, n - m, q, pow_p, carry);
00112   }
00113 public:
00114   template<typename C, typename M, typename I> static inline nat
00115   direct (I* c, nat n, const C& a, const M* pow_p) {
00116     I carry (0); nat o= direct (c, n, a, pow_p, carry); (void) carry;
00117     return o; }
00118 
00119 };
00120 
00121 
00122 
00123 
00124 
00125 template<typename C, typename S=std_base_dicho<C>,
00126          typename V=typename Base_dicho_variant(C) >
00127 struct base_dicho_transformer : public S {
00128 
00129   typedef typename S::modulus_base I;
00130   typedef modulus<I, typename S::modulus_base_variant> M;
00131   typedef implementation<base_transform,V> Base;
00132 
00133   vector<M> pow_p;
00134   M p;
00135 
00136 private:
00137   void extend (nat n) {
00138     while (((nat) 1) << N(pow_p) <= ((n + 1) >> 1))
00139       pow_p << M (square (* pow_p[N(pow_p) - 1])); }
00140 public:
00141 
00142   template<typename K> 
00143   inline base_dicho_transformer (const K& _p) : p(_p) {
00144     ASSERT (_p != 0, "invalid base");
00145     pow_p << p; }
00146 
00147   inline nat direct_transform (I* c, nat n, const C& a) {
00148     extend (n);
00149     return Base::direct (c, n, a, seg (pow_p)); }
00150 
00151   inline void inverse_transform (C& a, const I* c, nat n) {
00152     extend (n);
00153     Base::inverse (a, c, n, seg (pow_p)); }
00154 };
00155 
00156 } 
00157 #endif //__MMX__BASE_DICHO__HPP