00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_SPARSE_VECTOR_HPP
00014 #define __MMX_SPARSE_VECTOR_HPP
00015 #include <basix/pair.hpp>
00016 #include <basix/symbol.hpp>
00017
00019
00020 namespace mmx {
00021 struct eq_sparse_vector;
00022 #define TMPL_DEF \
00023 template<typename C, typename T, typename V= exact_eq_sparse_vector>
00024 #define TMPL template<typename C, typename T, typename V>
00025 #define Sparse_vector sparse_vector<C,T,V>
00026 #define Sparse_vector_rep sparse_vector_rep<C,T,V>
00027 #define Pair pair<T,C>
00028 TMPL class sparse_vector_rep;
00029 TMPL class sparse_vector;
00030 TMPL inline nat N (const Sparse_vector& v);
00031
00032
00033
00034
00035
00036 struct hard_eq_sparse_vector {
00037 typedef hard_eq_op key_op;
00038 typedef equal_op val_op;
00039 typedef hard_less_op l_op;
00040 };
00041
00042 struct exact_eq_sparse_vector {
00043 typedef exact_eq_op key_op;
00044 typedef equal_op val_op;
00045 typedef less_op l_op;
00046 };
00047
00048
00049
00050
00051
00052 TMPL_DEF
00053 class sparse_vector_rep REP_STRUCT {
00054 Pair* a;
00055 nat n;
00056 nat l;
00057
00058 public:
00059 inline sparse_vector_rep (Pair* a2, nat n2):
00060 a (a2), n (n2), l (n2) {}
00061 inline sparse_vector_rep (Pair* a2, nat n2, nat l2):
00062 a (a2), n (n2), l (l2) {}
00063 inline ~sparse_vector_rep () {
00064 mmx_delete<Pair > (a, l); }
00065
00066 friend class Sparse_vector;
00067 friend nat N LESSGTR (const Sparse_vector& v);
00068 };
00069
00070 TMPL_DEF
00071 class sparse_vector {
00072 INDIRECT_PROTO_3 (sparse_vector, sparse_vector_rep, C, T, V)
00073 inline sparse_vector () {
00074 rep= new Sparse_vector_rep (mmx_new<Pair > (0), 0); }
00075 inline sparse_vector (const T& k, const C& v) {
00076 typedef typename V::val_op Eq;
00077 static C zero; set_zero (zero);
00078 if (Eq::op (v, zero)) rep= new Sparse_vector_rep (mmx_new<Pair > (0), 0);
00079 else rep= new Sparse_vector_rep (mmx_new<Pair > (1, Pair (k, v)), 1); }
00080 inline sparse_vector (const Pair& p) {
00081 rep= new Sparse_vector_rep (mmx_new<Pair > (1, p), 1); }
00082 inline sparse_vector (Pair* a, nat n) {
00083 rep= new Sparse_vector_rep (a, n); }
00084 inline sparse_vector (Pair* a, nat n, nat l) {
00085 rep= new Sparse_vector_rep (a, n, l); }
00086 inline Pair operator [] (nat i) const {
00087 VERIFY (i<rep->n, "index out of range");
00088 return rep->a[i]; }
00089 inline Pair& operator [] (nat i) {
00090 VERIFY (i<rep->n, "index out of range");
00091 secure ();
00092 return rep->a[i]; }
00093 friend nat N LESSGTR (const Sparse_vector& v);
00094 };
00095 INDIRECT_IMPL_3 (sparse_vector, sparse_vector_rep,
00096 typename C, C, typename T, T, typename V, V)
00097
00098 TMPL inline nat N (const Sparse_vector& v) { return v.rep->n; }
00099
00100
00101
00102
00103
00104 TMPL syntactic
00105 flatten (const Sparse_vector& v) {
00106 vector<syntactic> r;
00107 for (nat i=0; i<N(v); i++)
00108 r << flatten (v[i]);
00109 return apply ("sparse_vector", r);
00110 }
00111
00112 TMPL
00113 struct binary_helper<Sparse_vector >:
00114 public void_binary_helper<Sparse_vector >
00115 {
00116 static inline string short_type_name () {
00117 return "Sv" * Short_type_name (C); }
00118 static inline generic full_type_name () {
00119 return gen ("Sparse_vector", Full_type_name (C)); }
00120 static inline nat size (const Sparse_vector& v) {
00121 return N(v); }
00122 static inline generic access (const Sparse_vector& v, nat i) {
00123 return as<generic> (v[i]); }
00124 static inline generic disassemble (const Sparse_vector& v) {
00125 vector<generic> r;
00126 for (nat i=0; i<N(v); i++) r << as<generic> (v[i]);
00127 return as<generic> (r); }
00128 static inline Sparse_vector assemble (const generic& x) {
00129 vector<generic> v= as<vector<generic> > (x);
00130 Pair* r= mmx_new<Pair > (N(v));
00131 for (nat i=0; i<N(v); i++) r[i]= as<Pair > (v);
00132 return Sparse_vector (r, N(v)); }
00133 static inline void write (const port& out, const Sparse_vector& v) {
00134 binary_write<nat> (out, N(v));
00135 for (nat i=0; i<N(v); i++)
00136 binary_write<Pair > (out, v[i]); }
00137 static Sparse_vector read (const port& in, nat n) {
00138 if (n == 0) return Sparse_vector ();
00139 else if (n == 1) return Sparse_vector (binary_read<Pair > (in));
00140 else return read (in, (n+1) >> 1) + read (in, n >> 1); }
00141 static inline Sparse_vector read (const port& in) {
00142 nat n= binary_read<nat> (in);
00143 return read (in, n); }
00144 };
00145
00146
00147
00148
00149
00150 template<typename Op, typename C, typename T, typename V> nat
00151 unary_hash (const Sparse_vector& v) {
00152 nat h=1236321;
00153 for (nat i=0; i<N(v); i++)
00154 h = (h << 5) ^ (h >> 27) ^ Op::op (v[i]);
00155 return h;
00156 }
00157
00158 template<typename Op, typename C, typename T, typename V> bool
00159 binary_test (const Sparse_vector& v1, const Sparse_vector& v2) {
00160 nat n1= N(v1), n2= N(v2);
00161 if (n1 != n2) return false;
00162 for (nat i=0; i<n1; i++)
00163 if (!Op::op (v1[i], v2[i]))
00164 return false;
00165 return true;
00166 }
00167
00168 TRUE_IDENTITY_OP_SUGAR(TMPL,Sparse_vector)
00169 EXACT_IDENTITY_OP_SUGAR(TMPL,Sparse_vector)
00170
00171
00172
00173
00174
00175 template<typename Op, typename C, typename T, typename V> Sparse_vector
00176 unary_map (const Sparse_vector& v) {
00177 typedef typename V::val_op Eq;
00178 static C zero; set_zero (zero);
00179 nat i, j, n= N(v);
00180 Pair* r= mmx_new<Pair > (n);
00181 for (i= j= 0; i < n; i++) {
00182 Pair e (v[i].x1, Op::op (v[i].x2));
00183 if (Eq::not_op (e.x2, zero)) r[j++]= e;
00184 }
00185 return Sparse_vector (r, j, n);
00186 }
00187
00188 TMPL inline Sparse_vector
00189 operator - (const Sparse_vector& t) {
00190 return unary_map<neg_op,C,T,V> (t); }
00191
00192
00193
00194
00195
00196 template<typename Op, typename C, typename T, typename V> Sparse_vector
00197 binary_map (const Sparse_vector& v1, const Sparse_vector& v2) {
00198 typedef typename V::val_op Eq;
00199 typedef typename V::l_op Less;
00200 static C zero; set_zero (zero);
00201 nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
00202 Pair* r= mmx_new<Pair > (n);
00203 for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
00204 if (Less::op (v1[i1].x1, v2[i2].x1)) {
00205 Pair e (v1[i1].x1, Op::op (v1[i1].x2, zero));
00206 i1++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00207 }
00208 else if (Less::op (v2[i2].x1, v1[i1].x1)) {
00209 Pair e (v2[i2].x1, Op::op (zero, v2[i2].x2));
00210 i2++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00211 }
00212 else {
00213 Pair e (v2[i2].x1, Op::op (v1[i1].x2, v2[i2].x2));
00214 i1++; i2++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00215 }
00216 }
00217 while (i1 < n1) {
00218 Pair e (v1[i1].x1, Op::op (v1[i1].x2, zero));
00219 i1++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00220 }
00221 while (i2 < n2) {
00222 Pair e (v2[i2].x1, Op::op (zero, v2[i2].x2));
00223 i2++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00224 }
00225 return Sparse_vector (r, j, n);
00226 }
00227
00228 template<typename Op, typename C, typename T, typename V> Sparse_vector
00229 binary_map_optimized (const Sparse_vector& v1, const Sparse_vector& v2) {
00230
00231 typedef typename V::val_op Eq;
00232 typedef typename V::l_op Less;
00233 static C zero; set_zero (zero);
00234 nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
00235 Pair* r= mmx_new<Pair > (n);
00236 for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
00237 if (Less::op (v1[i1].x1, v2[i2].x1)) {
00238 typedef typename Op::lop Lop;
00239 r[j++]= Pair (v1[i1].x1, Lop::op (v1[i1].x2)); i1++;
00240 }
00241 else if (Less::op (v2[i2].x1, v1[i1].x1)) {
00242 typedef typename Op::rop Rop;
00243 r[j++]= Pair (v2[i2].x1, Rop::op (v2[i2].x2)); i2++;
00244 }
00245 else {
00246 Pair e (v2[i2].x1, Op::op (v1[i1].x2, v2[i2].x2));
00247 i1++; i2++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00248 }
00249 }
00250 while (i1 < n1) {
00251 typedef typename Op::lop Lop;
00252 r[j++]= Pair (v1[i1].x1, Lop::op (v1[i1].x2)); i1++;
00253 }
00254 while (i2 < n2) {
00255 typedef typename Op::rop Rop;
00256 r[j++]= Pair (v2[i2].x1, Rop::op (v2[i2].x2)); i2++;
00257 }
00258 return Sparse_vector (r, j, n);
00259 }
00260
00261 TMPL inline Sparse_vector
00262 operator + (const Sparse_vector& t, const Sparse_vector& u) {
00263 return binary_map_optimized<add_op,C,T,V> (t, u); }
00264 TMPL inline Sparse_vector
00265 operator - (const Sparse_vector& t, const Sparse_vector& u) {
00266 return binary_map_optimized<sub_op,C,T,V> (t, u); }
00267 TMPL inline Sparse_vector
00268 operator * (const Sparse_vector& t, const Sparse_vector& u) {
00269 return binary_map<mul_op,C,T,V> (t, u); }
00270 TMPL inline Sparse_vector
00271 sup (const Sparse_vector& t, const Sparse_vector& u) {
00272 return binary_map<sup_op,C,T,V> (t, u); }
00273 TMPL inline Sparse_vector
00274 inf (const Sparse_vector& t, const Sparse_vector& u) {
00275 return binary_map<inf_op,C,T,V> (t, u); }
00276
00277
00278
00279
00280
00281 template<typename Op, typename C, typename T, typename V>
00282 Sparse_vector
00283 binary_map_scalar (const Sparse_vector& v, const C& x) {
00284 typedef typename V::val_op Eq;
00285 static C zero; set_zero (zero);
00286 nat i, j, n= N(v);
00287 Pair* r= mmx_new<Pair > (n);
00288 for (i=0, j=0; i<n; i++) {
00289 Pair e (v[i].x1, Op::op (v[i].x2, x));
00290 if (Eq::not_op (e.x2, zero)) r[j++]= e;
00291 }
00292 return Sparse_vector (r, j, n);
00293 }
00294
00295 TMPL inline Sparse_vector
00296 operator * (const Sparse_vector& t, const C& sc) {
00297 return binary_map_scalar<rmul_op> (t, sc); }
00298 TMPL inline Sparse_vector
00299 operator * (const C& sc, const Sparse_vector& t) {
00300 return binary_map_scalar<lmul_op> (t, sc); }
00301 TMPL inline Sparse_vector
00302 operator / (const Sparse_vector& t, const C& sc) {
00303 return binary_map_scalar<rdiv_op> (t, sc); }
00304 TMPL inline Sparse_vector
00305 operator / (const C& sc, const Sparse_vector& t) {
00306 return binary_map_scalar<ldiv_op> (t, sc); }
00307
00308
00309
00310
00311
00312 template<typename Mul, typename C, typename T, typename V> Sparse_vector
00313 _mul_add (const Sparse_vector& v1, const Sparse_vector& v2, const C& x) {
00314
00315 typedef typename V::val_op Eq;
00316 typedef typename V::l_op Less;
00317 static C zero; set_zero (zero);
00318 if (exact_eq (x, zero)) return v1;
00319 nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
00320 Pair* r= mmx_new<Pair > (n);
00321 for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
00322 if (Less::op (v1[i1].x1, v2[i2].x1)) {
00323 r[j++]= v1[i1]; i1++;
00324 }
00325 else if (Less::op (v2[i2].x1, v1[i1].x1)) {
00326 r[j++]= Pair (v2[i2].x1, Mul::op (x, v2[i2].x2)); i2++;
00327 }
00328 else {
00329 Pair e (v2[i2].x1, v1[i1].x2 + Mul::op (x, v2[i2].x2));
00330 i1++; i2++; if (Eq::not_op (e.x2, zero)) r[j++]= e;
00331 }
00332 }
00333 while (i1 < n1) {
00334 r[j++]= Pair (v1[i1].x1, v1[i1].x2); i1++;
00335 }
00336 while (i2 < n2) {
00337 r[j++]= Pair (v2[i2].x1, Mul::op (x, v2[i2].x2)); i2++;
00338 }
00339 return Sparse_vector (r, j, n);
00340 }
00341
00342 template<typename C, typename T, typename V> inline Sparse_vector
00343 sparse_mul_add (const Sparse_vector& v1, const Sparse_vector& v2, const C& x) {
00344 return _mul_add<mul_op> (v1, v2, x); }
00345
00346
00347
00348
00349
00350 #define TCM template<typename C, typename M>
00351
00352 TCM sparse_vector<C,symbol<M>,hard_eq_sparse_vector>
00353 make_formal_vector (const C& c, const M& m) {
00354 return sparse_vector<C,symbol<M>,hard_eq_sparse_vector> (symbol<M> (m), c);
00355 }
00356
00357 TCM C
00358 get_co (const sparse_vector<C,symbol<M>,hard_eq_sparse_vector>& v, nat i) {
00359 return v[i].x2;
00360 }
00361
00362 TCM M
00363 get_mo (const sparse_vector<C,symbol<M>,hard_eq_sparse_vector>& v, nat i) {
00364 return *v[i].x1;
00365 }
00366
00367 #undef TCM
00368
00369 #undef TMPL_DEF
00370 #undef TMPL
00371 #undef Sparse_vector
00372 #undef Sparse_vector_rep
00373 #undef Pair
00374 }
00375 #endif // __MMX_SPARSE_VECTOR_HPP
