00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__CRT_DICHO__HPP
00014 #define __MMX__CRT_DICHO__HPP
00015 #include <algebramix/crt_naive.hpp>
00016
00017 namespace mmx {
00018
00019
00020
00021
00022
00023 template<typename C>
00024 struct std_crt_dicho : std_crt_naive<C> {};
00025
00026
00027
00028
00029
00030 #define Crt_dicho_variant(C) crt_dicho_variant_helper<C>::CV
00031
00032 template<typename V>
00033 struct crt_dicho: public V {};
00034
00035 template<typename F, typename V, typename W>
00036 struct implementation<F,V,crt_dicho<W> >:
00037 public implementation<F,V,W> {};
00038
00039 template<typename C>
00040 struct crt_dicho_variant_helper {
00041 typedef crt_dicho<typename crt_naive_variant_helper<C>::CV> CV;
00042 };
00043
00044
00045
00046
00047
00048 template<typename V,typename W>
00049 struct implementation<crt_transform,V,crt_dicho<W> > :
00050 public implementation<crt_project,V> {
00051 typedef implementation<crt_project,V> Crt;
00052
00053 template<typename C, typename Modulus> static inline void
00054 down (C* c, const C& a, Modulus** tree, nat* k, nat depth, nat n) {
00055 if (n == 0) return;
00056 c[0]= Crt::mod (a, tree[depth-1][0]);
00057 if (depth == 1) return;
00058 for (nat i= depth-1; i > 0; i--) {
00059 if (k[i-1] & 1) c[k[i-1] - 1]= c[k[i] - 1];
00060 for (nat j= (k[i-1] >> 1) - 1; j + 1 > 0; j--) {
00061 c[(j << 1) + 1]= Crt::mod (c[j], tree[i-1][(j << 1) + 1]);
00062 c[(j << 1) ]= Crt::mod (c[j], tree[i-1][(j << 1) ]); } } }
00063
00064 template<typename C, typename Modulus> static inline void
00065 up (C& a, C* c, Modulus** tree, nat* k, nat depth, nat n) {
00066 if (n == 0) return;
00067 for (nat i= 0; i < depth - 1; i++) {
00068 for (nat j= 0; j < (k[i] >> 1); j++)
00069 c[j]= c[j << 1] * * tree[i][(j << 1) + 1] +
00070 c[(j << 1) + 1] * * tree[i][j << 1];
00071 if (k[i] & 1) c[k[i+1] - 1]= c[k[i] - 1]; }
00072 a= c[0]; }
00073
00074 template<typename C, typename J, typename Modulus> static inline void
00075 direct (J* c, const C& a,
00076 Modulus** tree, nat* k, nat depth, nat n, C* aux) {
00077 down (aux, a, tree, k, depth, n);
00078 for (nat i= 0; i < n; i++) c[i]= as<J> (aux[i]); }
00079
00080 template<typename C, typename Modulus> static inline void
00081 direct (C* c, const C& a,
00082 Modulus** tree, nat* k, nat depth, nat n, C* aux) {
00083 (void) aux; down (c, a, tree, k, depth, n); }
00084
00085 template<typename C, typename I, typename Modulus> static inline void
00086 combine (C& a, const I* c, Modulus** tree, nat* k, nat depth, nat n,
00087 C* aux) {
00088 if (n == 0) { a= C(0); return; }
00089 for (nat i= 0; i < n; i++) aux[i]= C (c[i]);
00090 up (a, aux, tree, k, depth, n); }
00091
00092 template<typename C, typename I, typename Modulus> static inline void
00093 inverse (C& a, const I* c,
00094 Modulus** tree, nat* k, nat depth,
00095 C* m, nat n, const Modulus& P, C* aux) {
00096 if (n == 0) { a= C(0); return; }
00097 for (nat i= 0; i < n; i++)
00098 aux[i]= Crt::mod (C (c[i]) * m[i], tree[0][i]);
00099 up (a, aux, tree, k, depth, n);
00100 a= Crt::mod (a, P); }
00101
00102 };
00103
00104
00105
00106
00107
00108 template<typename C, typename S=std_crt_dicho<C>,
00109 typename V=typename Crt_dicho_variant(C)>
00110 struct crt_dicho_transformer : public S {
00111 private:
00112 typedef typename S::modulus_base I;
00113 typedef modulus<I, typename S::modulus_base_variant> M;
00114 typedef modulus<C, typename S::modulus_variant> Modulus;
00115 typedef implementation<vector_linear, vector_naive> Vec;
00116 typedef implementation<crt_transform,V> Crt;
00117
00118 nat n;
00119 M* p;
00120 modulus<C>** tree;
00121 nat depth;
00122 nat* k;
00123 Modulus P;
00124 C H;
00125 C* m;
00126 C* aux;
00127
00128 inline void setup_direct (const M* _p, nat nn) {
00129 m= NULL;
00130 n= nn;
00131 if (n == 0) { depth= 0; P= 1; H= 0; return; }
00132 p= mmx_new<M> (n);
00133 for (nat i= 0; i < n; i++) p[i]= _p[i];
00134 depth= log_2 (next_power_of_two (n)) + 1;
00135 tree= mmx_new<Modulus*> (depth);
00136 k= mmx_new<nat> (depth);
00137 k[0]= n;
00138 tree[0]= mmx_new<Modulus> (k[0]);
00139 for (nat j= 0; j < k[0]; j++)
00140 tree[0][j]= C(* p[j]);
00141 for (nat i= 1; i < depth; i++) {
00142 k[i]= (k[i-1] >> 1) + (k[i-1] & 1);
00143 tree[i]= mmx_new<Modulus> (k[i]);
00144 for (nat j= 0; j < (k[i-1] >> 1); j++)
00145 tree[i][j]= * tree[i-1][j << 1] * * tree[i-1][(j << 1) + 1];
00146 if (k[i-1] & 1) tree[i][k[i] - 1]= * tree[i-1][k[i-1] - 1];
00147 }
00148 P= tree[depth-1][0];
00149 H= Crt::half (P); }
00150
00151 public:
00152
00153 inline void setup_inverse () {
00154 typedef implementation<vector_linear,vector_naive> Vec;
00155 if (n == 0 || m != NULL) return;
00156 m= mmx_new<C> (n);
00157 C c;
00158 Vec::set (m, C(1), n);
00159 Crt::up (c, m, tree, k, depth, n);
00160 Crt::down (m, c, tree, k, depth, n);
00161 for (nat i= 0; i < n; i++)
00162 inv_mod (m[i], tree[0][i]); }
00163
00164 inline crt_dicho_transformer (const vector<M>& p, bool lazy= true)
00165 : tree(NULL), k(NULL) {
00166 setup_direct (seg (p), N(p));
00167 aux= mmx_new<C> (n);
00168 if (! lazy) setup_inverse (); }
00169
00170 inline ~crt_dicho_transformer () {
00171 if (n == 0) return;
00172 mmx_delete<M> (p, n);
00173 mmx_delete<C> (aux, n);
00174 for (nat i= 0; i < depth; i++)
00175 mmx_delete<Modulus> (tree[i], k[i]);
00176 mmx_delete<Modulus*> (tree, depth);
00177 mmx_delete<nat> (k, depth);
00178 if (m != NULL) mmx_delete<C> (m, n); }
00179
00180 inline M operator[] (nat i) const {
00181 VERIFY (i < n, "index out of range");
00182 return p[i]; }
00183
00184 inline nat size () const {
00185 return n; }
00186
00187 inline Modulus product () const {
00188 return P; }
00189
00190 inline C comodulus (nat i) const {
00191 VERIFY (i < n, "index out of range");
00192 return *P / m[i]; }
00193
00194 inline void direct_transform (I* c, const C& a) const {
00195 C b= Crt::encode (a, P);
00196 Crt::direct (c, b, tree, k, depth, n, aux); }
00197
00198 inline void combine (C& a, const I* c) {
00199 Crt::combine (a, c, tree, k, depth, n, aux); }
00200
00201 inline void inverse_transform (C& a, const I* c) {
00202 setup_inverse ();
00203 Crt::inverse (a, c, tree, k, depth, m, n, P, aux);
00204 a= Crt::decode (a, P, H); }
00205 };
00206
00207 template<typename C, typename M, typename V> inline nat
00208 N (const crt_dicho_transformer<C,M,V>& crter) {
00209 return crter.size ();
00210 }
00211
00212 }
00213 #endif //__MMX__CRT_DICHO__HPP
