00001 #ifndef realroot_MPOLDSE_BERNSTEIN_C
00002 #define realroot_MPOLDSE_BERNSTEIN_C
00003
00004 namespace mmx {
00005
00006 namespace tensor
00007 {
00008 template< class C > inline void
00009 binoms( bernstein<C>& mpl ) {
00010 binoms( mpl.begin(), mpl.env, binomials<C>::default_ ) ;
00011 }
00012
00013 template < class C > inline void
00014 scale (bernstein < C > &mpl)
00015 {
00016 scale (mpl.begin (), mpl.env, binomials < C >::default_);
00017 };
00018
00019 template < class C > inline void
00020 uscale (bernstein < C > &mpl)
00021 {
00022 uscale (mpl.begin (), mpl.env, binomials < C >::default_);
00023 };
00024
00025 template < class C > inline void
00026 convertb2m ( monomials<C> & mpl ) {
00027 convertb2m (mpl.begin (), mpl.env, binomials < C >::default_);
00028 }
00029
00030 template < class C > inline void
00031 convertm2b (bernstein < C > &mpl) {
00032 convertm2b (mpl.begin (), mpl.env, binomials < C >::default_);
00033 };
00034
00035 template < class C , class V> inline void
00036 convertm2b (bernstein < C > &mpl, const V& bx) {
00037 convertm2b (mpl);
00038 for(unsigned i=0;i<bx.size();i+=2)
00039 restrict(mpl,i/2,bx[i],bx[i+1]);
00040 };
00041
00042 template < class C > void
00043 elevate (bernstein < C > &r, const eenv & elev)
00044 {
00045 bernstein < C > tmp (elev);
00046 tensor::binoms (tmp);
00047 scale (r);
00048 tensor::conv(r, tmp);
00049 uscale (r);
00050 };
00051
00052 template < class C > void
00053 add (bernstein < C > &mpl, const C & c)
00054 {
00055 vct::scadd (mpl.begin (), c, mpl.esz (), 1);
00056 };
00057
00058 template < class C > void
00059 add (bernstein < C > &r, const bernstein < C > &a, const C & c)
00060 {
00061 if (r.env != a.env)
00062 add (r = a, c);
00063 else
00064 vct::scadd (r.begin (), a.begin (), c, r.esz (), 1);
00065 };
00066
00067 template < class C > void
00068 add (bernstein < C > &r, const C & c, const bernstein < C > &a) {
00069 add(r,a,c);
00070 }
00071
00072 template<class C> void
00073 rewrite( bernstein<C>& r, const eenv& newe )
00074 {
00075 if ( newe == r.env ) return;
00076 eenv elenv = eenv::elevation (newe, r.env);
00077 elevate (r, elenv) ;
00078 };
00079
00080 template < class C > void
00081 add (bernstein < C > &r, const bernstein < C > &a)
00082 {
00083 eenv cm(eenv::common(r.env,a.env));
00084 if ( r.env != cm ) { rewrite(r,cm); };
00085 if ( a.env == cm )
00086 {
00087 vct::padd (r.begin (), a.begin (), r.esz (), 1, 1);
00088 }
00089 else
00090 {
00091 bernstein<C> tmp(a);
00092 rewrite(tmp,cm);
00093 vct::padd (r.begin (), tmp.begin (), r.esz (), 1, 1);
00094 };
00095 };
00096
00097
00098 template < class C > void
00099 add (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
00100 add (r=a, b);
00101 }
00102
00103
00104 template < class C > void
00105 sub (bernstein < C > &r, const bernstein < C > &a) {
00106 eenv cm(eenv::common(r.env,a.env));
00107 if ( r.env != cm ) { rewrite(r,cm); };
00108 if ( a.env == cm )
00109 {
00110 vct::psub (r.begin (), a.begin (), r.esz (), 1, 1);
00111 }
00112 else
00113 {
00114 bernstein<C> tmp(a);
00115 rewrite(tmp,cm);
00116 vct::psub (r.begin (), tmp.begin (), r.esz (), 1, 1);
00117 };
00118 }
00119
00120
00121 template < class C > void
00122 sub (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
00123 sub (r = a, b);
00124 }
00125
00126 template < class C > void
00127 sub (bernstein < C > &mpl, const C & c) {
00128 vct::scsub (mpl.begin (), c, mpl.esz (), 1);
00129 }
00130
00131 template < class C > void
00132 sub (bernstein < C > &mpl, const C & c, const bernstein < C > &a) {
00133 add(mpl,a,-c);
00134 }
00135
00136 template < class C > void
00137 sub (bernstein < C > &mpl, const bernstein < C > &a, const C& c) {
00138 add(mpl,a,-c);
00139 }
00140
00141 template < class C > void
00142 mul (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
00143 bernstein < C > a_ (a);
00144 scale (a_);
00145 bernstein < C > b_ (b);
00146 scale (b_);
00147 tensor::conv (r, a_, b_);
00148 uscale (r);
00149 };
00150
00151 template < class C > void
00152 mul (bernstein < C > &r, const bernstein < C > &a) {
00153 bernstein < C > a_;
00154 r.swap (a_);
00155 scale (a_);
00156 bernstein < C > b_ (a);
00157 scale (b_);
00158 tensor::conv (r, a_, b_);
00159 uscale (r);
00160 };
00161
00162 template < class C > inline void
00163 mul (bernstein < C > &r, const bernstein < C > &a, const C &c) {
00164 mul(r=a,c);
00165 };
00166
00167 template < class C > inline void
00168 mul (bernstein < C > &r, const C &c, const bernstein < C > &a) {
00169 mul(r=a,c);
00170 };
00171
00172 template < class C > void
00173 mul (bernstein < C > &r, const C &c) {
00174 vct::scmul(r.begin (),c, r.esz (), 1);
00175 }
00176
00177 template < class C > inline void
00178 div (bernstein < C > &r, const bernstein < C > &a, const C &c) {
00179 div(r=a,c);
00180 }
00181
00182 template < class C > bernstein < C >::bernstein (const eenv & e): base_type(e) {};
00183
00184 template < class C > bernstein < C >::bernstein (const eenv & e, const C & c):
00185 base_type(e,c){};
00186
00187 template < class C >
00188 bernstein < C >::bernstein (int nvr, const int *szs, const int *vrs): base_type(nvr,szs,vrs)
00189 {
00190 }
00191
00192 template < class C > bernstein < C >::bernstein (const C & x) : base_type(x) {}
00193
00194 template < class Result, class Coeff, class Parameters > inline
00195 void eval (Result & result, const bernstein < Coeff > & controls , const Parameters & parameters )
00196 {
00197 levalb (result, controls.begin (), controls.env,parameters);
00198 };
00199
00200 template < class Coeff> inline
00201 void eval (Coeff& result, const bernstein < Coeff > & controls , const Coeff & v)
00202 {
00203 Coeff tmp[controls.env.sz()];
00204 std::copy (controls.begin(), controls.end(), tmp);
00205 decasteljau (tmp , controls.env.sz(), v, 1);
00206 result=tmp[0];
00207 };
00208
00209 template < class Result, class Coeff, class Parameters > inline
00210 void eval (Result & result, const bernstein < Coeff > & controls , const Parameters & parameters , unsigned n) {
00211
00212 levalb( result, controls.begin (), controls.env, parameters );
00213 };
00214
00215 template < class Result, class Coeff, class Parameters > inline
00216 void subs0 (Result *result, const bernstein < Coeff > & controls , const Parameters & parameters )
00217 {
00218 subs0b( result, controls.begin (), controls.env, parameters );
00219 };
00220
00221 template < class C > void
00222 casteljau (bernstein < C > &a, bernstein < C > &b, const C & t, int v)
00223 {
00224 if (a.env != b.env)
00225 realloc (a, b.env);
00226 casteljau (a.begin (), b.begin (), t, v);
00227 };
00228
00229 template < class C > void
00230 casteljau (bernstein < C > &left, bernstein < C > &right,
00231 const bernstein < C > &source, const C & t, int v)
00232 {
00233 casteljau (left, right = source, t, v);
00234 };
00235
00236
00237 template < class C > void
00238 casteljau (bernstein < C > &a, bernstein < C > &b, int v)
00239 {
00240 if (a.env != b.env)
00241 realloc (a, b.env);
00242 casteljau (a.begin (), b.begin (), v);
00243 };
00244
00245 template < class C > void
00246 casteljau (bernstein < C > &left, bernstein < C > &right,
00247 const bernstein < C > &source, int v)
00248 {
00249 casteljau (left, right = source, v);
00250 };
00251
00252 template < class C >
00253 std::ostream & operator<< (std::ostream & o, const bernstein < C > &mpl)
00254 {
00255 o << "bernstein object:\n";
00256 o << "\t(addresses) o: " << &mpl << " e: " << mpl.env.data << "d: " << mpl.
00257 begin () << std::endl;
00258 o << "\t(eenv)" << mpl.env << std::endl;
00259 o << "\t(data)";
00260 array::print (o, mpl.data);
00261 o << std::endl;
00262 o << ":object bernstein\n";
00263 return o;
00264 };
00265
00266 template<class C> void
00267 assign( monomials<C>& monoms, const bernstein<C>& controls )
00268 {
00269 monoms = controls;
00270 tensor::convertb2m(monoms);
00271 };
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282 template < class C > void
00283 cfdump (std::ostream & o, const bernstein < C > &mpl)
00284 {
00285 for (const C * a = mpl.begin (); a != mpl.end (); o << " " << *a++){}
00286 };
00287
00288 template <class OSTREAM, class C > void
00289 mprint (OSTREAM & o, const bernstein < C > &mpl)
00290 {
00291 monomials < C > mp (0);
00292 convert (mp, mpl);
00293 o << mp;
00294
00295 };
00296
00297 template<class OSTREAM, class C > void
00298 print(OSTREAM& o , const bernstein<C>& mpl ) {
00299 monomials<C> mp(mpl);
00300 convertb2m(mp);
00301 print(o,mp);
00302 };
00303
00304 template<class OSTREAM, class C > void inline
00305 print(OSTREAM& o , const bernstein<C>& mpl , const variables& v) {
00306 monomials<C> mp(mpl);
00307 convertb2m(mp);
00308 print(o,mp,v);
00309 };
00310
00311 template < class C > bool bernstein < C >::operator== (const bernstein & mpl) const
00312 {
00313 return base_type::operator==(mpl);
00314 };
00315
00316 template<class C> inline
00317 void diff( bernstein<C>& res, const bernstein<C>& src, int v )
00318 {
00319 clear( res, eenv::diff(src.env,v) );
00320 bdiff( res.begin(), src.env, src.begin(), v );
00321 };
00322
00323 template<class C>
00324 bernstein<C>& bernstein<C>::operator=( const C& x )
00325 {
00326 this->~bernstein<C>();
00327 new (this) bernstein<C>(x);
00328 return *this;
00329 };
00330
00331 template<class C,class T>
00332 void split( bernstein<C>& l, bernstein<C>& r, int v, const T& t )
00333 {
00334 if ( l.env != r.env ) { realloc(l,r.env); };
00335 bsplit( l.begin(), r.begin(), r.env, t, r.env.localv(v) );
00336 };
00337
00338
00339 template<class C> inline
00340 void split( bernstein<C>& l, bernstein<C>& r, int v )
00341 {
00342 if ( l.env != r.env ) { realloc(l,r.env); };
00343 bsplit2( l.begin(), r.begin(), r.env, r.env.localv(v) );
00344 };
00345
00346 template<class C,class T>
00347 void restrict( bernstein<C>& rst, int v, const T& a, const T& b ) {
00348 v = rst.env.vridx(v);
00349 if (v < rst.env.nvr())
00350 brestrict( rst.begin(), rst.env, v, a, b );
00351 };
00352
00353 template<class C> inline
00354 void lface ( bernstein<C>& f, const bernstein<C>& src, int v )
00355 {
00356 lface( f.begin(), src.env, src.begin(), v );
00357 };
00358
00359
00360 template<class C> inline
00361 void rface ( bernstein<C>& f, const bernstein<C>& src, int v )
00362 {
00363 rface( f.begin(), src.env, src.begin(), v );
00364 };
00365
00366
00367 template<class C> inline
00368 void face ( bernstein<C>& f, const bernstein<C>& src, int v, int n )
00369 {
00370 assert(n==0||n==1);
00371 eenv fenv(face_env(src.env,src.env.localv(v)));
00372 if ( f.env != fenv ) realloc( f, fenv );
00373 if ( n ) rface( f, src, src.env.localv(v) );
00374 else lface( f, src, src.env.localv(v) );
00375 };
00376
00377
00378 }
00379
00380
00381 namespace let {
00382 template<class C, class U> void assign(C& p, const U& q);
00383
00384 template<class C, class U>
00385 void assign(tensor::bernstein<C>& p, const tensor::bernstein<U>& q)
00386 {
00387 p.env=q.env;
00388 p.data = typename tensor::bernstein<C>::vector_type(q.data.size());
00389 for(unsigned i=0;i<q.data.size();i++) assign(p.data[i],q.data[i]);
00390 }
00391
00392 template<class C, class D> inline void
00393 assign( tensor::bernstein<C>& b, const tensor::monomials<D>& m)
00394 {
00395 b = tensor::bernstein<C>(m.env);
00396 for(int i=0;i<m.size();i++) b[i] =as<C>(m[i]);
00397 tensor::convertm2b(b);
00398 }
00399
00400 template<class C, class D, class O> inline void
00401 assign( tensor::bernstein<C>& b, const sparse::monomial_seq<D,O>& m)
00402 {
00403 b = tensor::bernstein<C>(m);
00404 }
00405 template<class C, class D, class O, class DOM> inline void
00406 assign( tensor::bernstein<C>& b, const sparse::monomial_seq<D,O>& m, const DOM& dom)
00407 {
00408 b = tensor::bernstein<C>(m,dom);
00409 }
00410
00411 }
00412
00413 }
00414
00415 #include "realroot/tensor_convert.hpp"
00416
00417 #endif
