00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN
00019 #define DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN
00020
00022
00031 template<class BASIS>
00032 class algebra : public sparse_vector<BASIS>
00033 {
00034 public:
00036 typedef sparse_vector<BASIS> VECT;
00038 typedef typename VECT::iterator iterator;
00040 typedef typename VECT::const_iterator const_iterator;
00042 typedef typename VECT::KEY KEY;
00044 typedef typename VECT::SCALAR SCALAR;
00046 typedef typename VECT::RATIONAL RATIONAL;
00047
00048 static const DEG MAX_DEGREE = BASIS::MAX_DEGREE;
00049
00050 template <class Transform, unsigned DEPTH1>
00051 inline void triangularbufferedmultiplyandcombine(const algebra& rhs, algebra& result, Transform fn) const
00052 {
00053
00054 std::vector<std::pair<KEY, SCALAR> > buffer(rhs.begin(), rhs.end());
00055 const_iterator i;
00056
00057
00058 std::vector<std::vector<std::pair<KEY, SCALAR> >:: const_iterator>
00059 iterators(DEPTH1 + 1,buffer.end());
00060 unsigned deg = 0;
00061 for (std::vector<std::pair<KEY, SCALAR> >:: const_iterator j0 = buffer.begin(); j0 != buffer.end(); j0++)
00062 {
00063 DEG d = basis.degree(j0->first);
00064 assert(d >= deg && d <= DEPTH1);
00065 while (deg < d)
00066 iterators[deg++] = j0;
00067
00068 }
00069
00070 std::vector<std::pair<KEY, SCALAR> >:: const_iterator j, jEnd;
00071 for (i = VECT::begin(); i != VECT::end(); ++i)
00072 {
00073 const KEY& k = i->first;
00074 unsigned rhdegree = DEPTH1 - basis.degree(k);
00075 std::vector<std::pair<KEY, SCALAR> >:: const_iterator&
00076 jEnd = iterators[rhdegree];
00077 for (j = buffer.begin(); j != jEnd; ++j)
00078 result.add_scal_prod(VECT::basis.prod(i->first, j->first),
00079 fn(i->second * j->second));
00080 }
00081
00082 }
00083
00084 template <class Transform>
00085 inline void squarebufferedmultiplyandcombine(const algebra& rhs, algebra& result, Transform fn) const
00086 {
00087
00088 std::vector<std::pair<KEY, SCALAR> > buffer(rhs.begin(), rhs.end());
00089 const_iterator i;
00090
00091
00092 std::vector<std::pair<KEY, SCALAR> >:: const_iterator j;
00093 for (i = VECT::begin(); i != VECT::end(); ++i)
00094 {
00095 for (j = buffer.begin(); j != buffer.end(); ++j)
00096 result.add_scal_prod(VECT::basis.prod(i->first, j->first),
00097 fn(i->second * j->second));
00098 }
00099
00100 }
00101
00102
00103
00105 struct scalar_passthrough
00106 {
00107 SCALAR operator()(const SCALAR& arg)
00108 {
00109 return arg;
00110 }
00111 };
00113 struct scalar_minus
00114 {
00115 SCALAR operator()(const SCALAR& arg)
00116 {
00117 return -arg;
00118 }
00119 };
00121 struct scalar_post_mult
00122 {
00123 private:
00124 SCALAR mFactor;
00125 public:
00126 scalar_post_mult(const SCALAR Factor = one)
00127 : mFactor(Factor) {}
00128 SCALAR operator()(const SCALAR arg)
00129 {
00130 return arg * mFactor;
00131 }
00132 };
00134 struct rational_post_div
00135 {
00136 private:
00137 SCALAR mFactor;
00138 public:
00139 rational_post_div(const RATIONAL Factor = one)
00140 : mFactor(one / Factor) {}
00141 SCALAR operator()(const SCALAR arg)
00142 {
00143 return arg * mFactor;
00144 }
00145 };
00147 template <unsigned DEPTH1>
00148 inline void bufferedmultiplyandadd(const algebra& rhs, algebra& result) const
00149 {
00150 scalar_passthrough fn;
00151 triangularbufferedmultiplyandcombine<scalar_passthrough, DEPTH1>(rhs, result, fn);
00152 }
00154 template <>
00155 inline void bufferedmultiplyandadd<0>(const algebra& rhs, algebra& result) const
00156 {
00157 scalar_passthrough fn;
00158 squarebufferedmultiplyandcombine(rhs, result, fn);
00159 }
00161 template <unsigned DEPTH1>
00162 inline void bufferedmultiplyandsub(const algebra& rhs, algebra& result) const
00163 {
00164 scalar_minus fn;
00165 triangularbufferedmultiplyandcombine<scalar_minus, DEPTH1>(rhs, result, fn);
00166 }
00167
00169 template <>
00170 inline void bufferedmultiplyandsub<0>(const algebra& rhs, algebra& result) const
00171 {
00172 scalar_minus fn;
00173 squarebufferedmultiplyandcombine(rhs, result, fn);
00174 }
00175
00176 struct wrapscalar
00177 {
00178 const SCALAR& hidden;
00179 wrapscalar(const SCALAR& s)
00180 : hidden(s) {}
00181 };
00182 struct wraprational
00183 {
00184 const RATIONAL& hidden;
00185 wraprational(const RATIONAL& s)
00186 : hidden(s) {}
00187 };
00189 template <unsigned DEPTH1>
00190 inline void bufferedmultiplyandsmult(const wrapscalar& ss, const algebra& rhs, algebra& result) const
00191 {
00192 scalar_post_mult fn(ss.hidden);
00193 triangularbufferedmultiplyandcombine<scalar_post_mult, DEPTH1>(rhs, result, fn);
00194 }
00195
00197 template <>
00198 inline void bufferedmultiplyandsmult<0>(const wrapscalar& ss, const algebra& rhs, algebra& result) const
00199 {
00200 scalar_post_mult fn(ss.hidden);
00201 squarebufferedmultiplyandcombine(rhs, result, fn);
00202 }
00203
00205 template <unsigned DEPTH1>
00206 inline void bufferedmultiplyandsdiv(const algebra& rhs, const wraprational& ss, algebra& result) const
00207 {
00208 rational_post_div fn(ss.hidden);
00209 triangularbufferedmultiplyandcombine<rational_post_div, DEPTH1>(rhs, result, fn);
00210 }
00211
00213 template <>
00214 inline void bufferedmultiplyandsdiv<0>(const algebra& rhs, const wraprational& ss, algebra& result) const
00215 {
00216 rational_post_div fn(ss.hidden);
00217 squarebufferedmultiplyandcombine(rhs, result, fn);
00218 }
00219
00220 public:
00222
00225 algebra(void) {}
00227 algebra(const algebra& a)
00228 : VECT(a) {}
00230 algebra(const VECT& v)
00231 : VECT(v) {}
00233 explicit algebra(const KEY& k, const SCALAR& s = VECT::one)
00234 : VECT(k, s) {}
00235 public:
00237 inline algebra& operator*=(const SCALAR& s)
00238 {
00239 VECT::operator *= (s);
00240 return *this;
00241 }
00243 inline algebra& operator/=(const RATIONAL& s)
00244 {
00245 VECT::operator /= (s);
00246 return *this;
00247 }
00249 inline __DECLARE_BINARY_OPERATOR(algebra,*,*=,SCALAR)
00251 inline __DECLARE_BINARY_OPERATOR(algebra,/,/=,SCALAR)
00253 inline __DECLARE_BINARY_OPERATOR(algebra,+,+=,algebra)
00255 inline __DECLARE_BINARY_OPERATOR(algebra,-,-=,algebra)
00257 inline __DECLARE_UNARY_OPERATOR(algebra,-,-,VECT);
00259 inline algebra& operator*=(const algebra& rhs)
00260 {
00261 algebra result;
00262 bufferedmultiplyandadd<MAX_DEGREE>(rhs, result);
00263 swap(result);
00264 return *this;
00265 }
00267 inline __DECLARE_BINARY_OPERATOR(algebra,*,*=,algebra);
00269 inline algebra& add_mul(const algebra& a, const algebra& b)
00270 {
00271 a.bufferedmultiplyandadd<MAX_DEGREE>(b, *this);
00272 return *this;
00273 }
00275 inline algebra& sub_mul(const algebra& a, const algebra& b)
00276 {
00277 a.bufferedmultiplyandsub<MAX_DEGREE>(b, *this);
00278 return *this;
00279 }
00281 inline algebra& mul_scal_prod(const algebra& rhs, const SCALAR& s)
00282 {
00283 algebra result;
00284 bufferedmultiplyandsmult(rhs, wrapscalar(s), result);
00285 swap(result);
00286 return *this;
00287 }
00289 inline algebra& mul_scal_div(const algebra& rhs, const RATIONAL& s)
00290 {
00291 algebra result;
00292 bufferedmultiplyandsdiv<MAX_DEGREE>(rhs, wraprational(s), result);
00293 swap(result);
00294 return *this;
00295 }
00297 inline friend algebra commutator(const algebra& a, const algebra& b)
00298 {
00299 algebra result;
00300 a.bufferedmultiplyandadd<MAX_DEGREE>(b, result);
00301 b.bufferedmultiplyandsub<MAX_DEGREE>(a, result);
00302 return result;
00303 }
00305 inline algebra truncate(const DEG min, const DEG max) const
00306 {
00307 algebra result;
00308 const_iterator i;
00309 for (i = VECT::begin(); i != VECT::end(); ++i)
00310 if ((VECT::basis.degree(i->first) >= min) && (VECT::basis.degree(i->first) <= max))
00311 result[i->first] = i->second;
00312 return result;
00313 }
00315 inline DEG degree(void) const
00316 {
00317 DEG result(0);
00318 const_iterator i;
00319 for (i = VECT::begin(); i != VECT::end(); ++i)
00320 result = std::max(result, VECT::basis.degree(i->first));
00321 return result;
00322 }
00323 };
00324
00325
00326 #endif // DJC_COROPA_LIBALGEBRA_ALGEBRAH_SEEN
00327
00328