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