00001 #ifndef realroot_MPOLDSE_VCTOPS_H
00002 #define realroot_MPOLDSE_VCTOPS_H
00003 #include <iostream>
00004 #include <realroot/texp_sup.hpp>
00005 #include <realroot/texp_ringof.hpp>
00006
00007 namespace mmx {
00008
00009 namespace vct
00010 {
00011
00012 template<typename A, typename B> inline
00013 void padd( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1 )
00014 {
00015 for ( A * ea = a; ea != a+sz*sta; *ea += *b , ea += sta, b += stb ) {}
00016 };
00017
00018 template<typename A, typename B, typename C> inline
00019 void padd( A * a, const B * b, const C * c, unsigned sz,
00020 int sta = 1, int stb = 1, int stc = 1 )
00021 {
00022 for ( A * ea = a; ea != a+sz*sta; *ea = *b+*c, ea += sta, b += stb, c += stc ) ;
00023 };
00024
00025 template<typename A, typename B> inline
00026 void psub( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1 )
00027 {
00028 for ( A * ea = a; ea != a+sz*sta; *ea -= *b, ea += sta, b += stb ) ;
00029 };
00030
00031 template<typename A, typename B, typename C> inline
00032 void psub( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1 )
00033 {
00034 for ( A * ea = a; ea != a+sz*sta; *ea = *b-*c, ea += sta, b += stb, c += stc ) ;
00035 };
00036
00037 template<typename A, typename B> inline
00038 void pmul( A * a, const B * b, int sz, int sta = 1, int stb = 1 )
00039 {
00040 for ( A * ea = a; ea != a+sz*sta; *ea *= *b, ea += sta, b += stb ) ;
00041 };
00042
00043 template<typename A, typename B, typename C> inline
00044 void pmul( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1 )
00045 {
00046 for ( A * ea = a; ea != a+sz*sta; *ea = *b**c, ea += sta, b += stb, c += stc ) ;
00047 };
00048
00049 template<typename A, typename B> inline
00050 void pdiv( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1 )
00051 {
00052 for ( A * ea = a; ea != a+sz*sta; *ea /= *b, ea += sta, b += stb ) ;
00053 };
00054
00055 template<typename A, typename B, typename C> inline
00056 void pdiv( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1 )
00057 {
00058 for ( A * ea = a; ea != a+sz*sta; *ea = *b/(*c), ea += sta, b += stb, c += stc ) ;
00059 };
00060
00061 template<typename A, typename B> inline
00062 void scadd( A * a, const B& b, unsigned n, int s = 1 )
00063 {
00064 for ( A * ea = a; ea != a + n*s; *a += b, ea += s ) ;
00065 };
00066
00067 template<typename A, typename B, typename C> inline
00068 void scadd( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
00069 {
00070 for ( A * ea = a; ea != a + n*sa; *ea = *b + c, ea += sa, b += sb ) ;
00071 };
00072
00073 template<typename A, typename B> inline
00074 void scsub( A * a, const B& b, unsigned n, int s = 1 )
00075 {
00076 for ( A * ea = a; ea != a + n*s; *ea -= b, ea += s ) ;
00077 };
00078
00079 template<typename A, typename B> inline
00080 void scsub( const B& b, A * a, unsigned n, int s = 1 )
00081 {
00082 for ( A * ea = a; ea != a+n*s; *ea = b-*ea, ea += s ) ;
00083 }
00084
00085 template<typename A, typename B, typename C> inline
00086 void scsub( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
00087 {
00088 for ( A * ea = a; ea != a+n*sa; *ea = *b-c, ea += sa, b += sb ) ;
00089 }
00090
00091 template<typename A, typename B, typename C> inline
00092 void scsub( A * a, const B& b, const C * c, unsigned n, int sa = 1, int sb = 1 )
00093 {
00094 for ( A * ea = a; ea != a+n*sa; *ea = c-*b, ea += sa, b += sb ) ;
00095 };
00096
00097 template<typename A, typename B> inline
00098 void scmul( A * a, const B& b, unsigned n, int s = 1 )
00099 {
00100 for ( A * ea = a; ea != a+n*s; *ea *= b, ea += s ) ;
00101 }
00102
00103 template<typename A, typename B, typename C> inline
00104 void scmul( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
00105 {
00106 for ( A * ea = a; ea != a + n*sa; *ea = *b*c, ea += sa, b += sb ) ;
00107 }
00108
00109 template<typename A, typename B, typename C> inline
00110 void scmul( A * a, const A * _a, const B * b, const C& c )
00111 {
00112 for ( A * ea = a; ea != _a; *ea++ = *b++*c ) ;
00113 }
00114
00115 template<typename A, typename B> inline
00116 void scdiv( A * a, const B& b, unsigned n, int s = 1 )
00117 {
00118 for ( A * ea = a; ea != a+n*s; *ea /= b, ea += s ) ;
00119 }
00120
00121 template<typename A, typename B> inline
00122 void scdiv( const B& b, A * a, unsigned n, int s = 1 )
00123 {
00124 for ( A * ea = a; ea != a+n*s; *ea = b/(*ea), ea += s ) ;
00125 };
00126
00127 template<typename A, typename B, typename C> inline
00128 void scdiv( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
00129 {
00130 for ( A * ea = a; ea != a + n*sa; *ea = *b/c, ea += sa, b += sb ) ;
00131 }
00132
00133 template<typename A, typename B, typename C> inline
00134 void scdiv( A * a, const C& c, const B * b , unsigned n, int sa = 1, int sb = 1 )
00135 {
00136 for ( A * ea = a; ea != a + n*sa; *ea = c/(*b), ea += sa, b += sb ) ;
00137 }
00138
00139 template<typename A, typename B> inline
00140 void copy ( A * a, const B * b, unsigned n )
00141 {
00142 std::copy(b,b+n,a) ;
00143 };
00144
00145 template<typename A, typename B> inline
00146 void copy ( A * a, const B * b, unsigned n, int sa, int sb )
00147 {
00148 for ( A * ea = a; ea != a + n*sa; *ea = *b, ea += sa, b += sb ) ;
00149 };
00150
00151 template<class C >
00152 void fill( C * a, unsigned n, int s, const C& x )
00153 {
00154 for ( C * pa = a; pa != a+n*s; *pa = x, pa += s ) ;
00155 };
00156
00157 template<typename A, typename B> inline
00158 void icopy( A * a, unsigned * aadd, unsigned nadd, const B * b )
00159 {
00160 for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] = *b++ ) ;
00161 };
00162
00163 template<typename A, typename B> inline
00164 void ipadd( A * a, unsigned * aadd, unsigned nadd, const B * b )
00165 {
00166 for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] += *b++ ) ;
00167 };
00168
00169 template<typename A, typename B> inline
00170 void ipsub( A * a, unsigned * aadd, unsigned nadd, const B * b )
00171 {
00172 for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] -= *b++ ) ;
00173 };
00174
00175
00176 template<typename A, typename B> inline
00177 void accmin( A& mn, B const * v, unsigned n, int s = 1 )
00178 {
00179 const B * p;
00180 for ( p = v; p != v + n*s; p += s )
00181 if ( mn > *p ) mn = *p;
00182 };
00183
00184 template<typename A, typename B> inline
00185 void min( A& mn, B const * v, unsigned n, int s = 1 )
00186 {
00187 mn = v[0];
00188 accmin(mn,v+s,n-1,s) ;
00189 };
00190
00191 template<typename A, typename B> inline
00192 void accmax( A& mx, const B * v, unsigned n, int s = 1 )
00193 {
00194 const B * p;
00195 for ( p = v; p != v+n*s; p += s )
00196 if ( mx < *p ) mx = *p;
00197 };
00198
00199 template<typename A, typename B> inline
00200 void max( A& mx, B const * v, unsigned n, int s = 1 )
00201 {
00202 mx = v[0];
00203 accmax(mx,v,n,s) ;
00204 };
00205 template<typename A> inline
00206 void print( const A * p, unsigned n, int st = 1,
00207 std::ostream& o = std::cout )
00208 {
00209 const A * i;
00210 o << "[ ";
00211 for ( i = p; i != p + (n-1)*st; i += st )
00212 o << *i << ", ";
00213 o << *i << " ]";
00214 };
00215
00228 template<class C>
00229 void conv( C * r,
00230 unsigned nr, int str,
00231 const
00232 C * a,
00233 unsigned na, int sta,
00234 const
00235 C * b,
00236 unsigned nb, int stb )
00237 {
00238 if ( na < nb ) { conv(r,nr,str,b,nb,stb,a,na,sta); return; };
00239 int k, l;
00240 C * er;
00241 const C * ea, * eb, * eeb, * eea;
00242 er = r;
00243 eb = b;
00244
00245 for ( er = r, eb = b, k = 0; k < std::min(nb,nr); k ++, er += str, eb += stb )
00246 for ( eeb = eb, ea = a, *er = 0, l = 0;
00247 l <= k;
00248 *er += *ea**eeb, l ++, eeb -= stb, ea += sta ) ;
00249
00250 for ( ; k < std::min(na,nr) ; k ++, er += str )
00251 for ( *er = 0,
00252 l = k-nb+1, eeb = b +(nb-1)*stb , ea = a + l*sta;
00253 l <= k;
00254 *er+= *ea**eeb,
00255 l ++, eeb -= stb, ea += sta ) ;
00256
00257 for ( ; k < nr; k ++, er += str )
00258 for ( *er = 0,
00259 l = k-nb+1, eeb = b+(nb-1)*stb, ea = a + l*sta;
00260 l < na;
00261 *er += *ea**eeb,
00262 l ++, eeb -= stb, ea += sta ) ;
00263 };
00264
00276 template<class C>
00277 void accconv( C * r, int str,
00278 const C * a, unsigned na, int sta,
00279 const C * b, unsigned nb, int stb )
00280 {
00281 C * er;
00282 const C * ea, * eb;
00283 for ( ea = a; ea != a + na*sta; ea += sta, r += str )
00284 for ( er = r, eb = b; eb != b + nb*stb; eb += stb, er += str )
00285 *er += *ea**eb;
00286 };
00287
00299 template<class C>
00300 void conv( C * r, int str,
00301 const C * a, unsigned na, int sta,
00302 const C * b, unsigned nb, int stb )
00303 {
00304 C * er;
00305 for ( er = r; er != r + (na+nb-1)*str; *er++ = 0 ) ;
00306 accconv(r,str,a,na,sta,b,nb,stb);
00307 };
00308
00316 template<class C> inline
00317 void minmax( C& m, C& M, const C * a, unsigned n, unsigned s = 1 )
00318 {
00319 for ( const C * e = a + n*s-s, m = M = *e; a < e; a += 2*s )
00320 {
00321 if ( *a > *(a+s) )
00322 {
00323 if ( *a > M ) M = *a;
00324 if ( *(a+s) < m ) m = *(a+s) ;
00325 }
00326 else
00327 {
00328 if ( *a > M ) M = *a;
00329 if ( *(a+s) < m ) m = *(a+s) ;
00330 };
00331 };
00332 };
00333
00334 template<class C> inline
00335 unsigned mmerge( C * r, const C* a, const C* b,
00336 unsigned na, unsigned nb )
00337 {
00338 C * er(r);
00339 const C * ea(a), * eb(b);
00340 a += na;
00341 b += nb;
00342 for (;;)
00343 {
00344 while( ea < a && *ea < *eb ) *er++ = *ea++;
00345 if ( ea == a ) { std::copy(eb,b,er) ; break; };
00346 if ( *ea == *eb ) { *er++ = *ea++; eb++;};
00347 while( eb < b && *eb < *ea ) *er++ = *eb++;
00348 if ( eb == b ) { std::copy(ea,a,er); break; };
00349 if ( *eb == *ea ) { *er++ = *eb++; ea++;};
00350 };
00351 return er-r;
00352 };
00353
00354 template<class C> inline
00355 void mmerge( C * r, int str,
00356 const C * a, unsigned na, int sta,
00357 const C * b, unsigned nb, int stb )
00358 {
00359 C * er(r);
00360 const C * ea(a), * eb(b);
00361
00362 a += na*sta;
00363 b += nb*stb;
00364
00365 for (;;)
00366 {
00367 for (;ea != a && *ea < *eb; *er = *ea, er += str, ea += sta ) ;
00368
00369 if ( ea == a )
00370 {
00371 for ( ; eb != b; *er = *eb, er += str, eb += stb ) ;
00372 break;
00373 };
00374
00375 if ( *ea == *eb )
00376 {
00377 *er = *ea;
00378 er += str;
00379 ea += sta;
00380 eb += stb;
00381 };
00382
00383 for (;eb != b && *eb < *ea; *er = *eb, er += str, eb += stb ) ;
00384
00385 if ( eb == b )
00386 {
00387 for ( ; ea != a; *er = *ea, er += str, ea += sta ) ;
00388 break;
00389 };
00390
00391 if ( *eb == *ea )
00392 {
00393 *er = *eb;
00394 er += str;
00395 eb += stb;
00396 ea += sta;
00397 };
00398 };
00399 return (er - r)/str;
00400 };
00401
00402
00403
00404 template < class C > static
00405 void shift (C * p, const C & c, unsigned n, int is = 1)
00406 {
00407 int j, k, s;
00408 for (s = n, j = 0; j <= s - 2; j++)
00409 for (k = s - 2; k >= j; p[k * is] += c * p[k * is + is], k--) ;
00410 };
00411
00412 template < class C > static
00413 void scale (C * p, const C & s, unsigned n, int is = 1)
00414 {
00415 C pw = s;
00416 for (C * cp = p + is; cp != p + n * is; *cp *= pw, pw *= s, cp += is) ;
00417 };
00418
00419 template<class O, class C, class I> inline
00420 void horner( O& o, C const * const c, unsigned n, const I& i, int s = 1 )
00421 {
00422 o = O(0);
00423 for ( int p = (n-1)*s; p != 0; o += as<O>(c[p]), o *= i, p -= s ) ;
00424 o+= as<O>(c[0]);
00425 };
00426
00427 template<typename C, typename I> inline typename texp::sup<C,I>::T
00428 horner( C const * const c, unsigned n, const I& i, int s = 1 )
00429 {
00430 typename texp::ringof<C,I>::T o(0);
00431 horner(o,c,n,i,s);
00432 return o;
00433 };
00434
00435 template<class O, class C, class I> inline
00436 void hhorner( O& o, C const * const c, unsigned n, const I& i, const I&i0, int s = 1 )
00437 {
00438 o = O(c[(n-1)*s]);
00439 for ( int p = (n-2)*s; p >= 0; o *= i, o += O(c[p])*i0, p -= s ) ;
00440 };
00441
00442
00443 template<typename C, typename I> inline typename texp::sup<C,I>::T
00444 hhorner( C const * const c, unsigned n, const I& i, const I&i0, int s = 1 )
00445 {
00446 typename texp::ringof<C,I>::T o(0);
00447 hhorner(o,c,n,i,i0,s) ;
00448 return o;
00449 };
00450
00451
00452 template<typename C> inline
00453 void diff( C * dst, C const * const src, unsigned sz, int st = 1 )
00454 {
00455 for ( unsigned i = 0; i < sz-1; dst[i*st] = (i+1)*src[(i+1)*st], i++ ) ;
00456 };
00457
00458 template<typename C, typename I>
00459 void decasteljau( C * c, unsigned n, const I& i, int s = 1 )
00460 {
00461 C *ec, *p;
00462 for ( ec = c + (n-1)*s; ec != c; ec -= s )
00463 for ( p = c; p != ec; p += s )
00464 *p = (1.0-i)**p+i**(p+s) ;
00465 };
00466
00467 template<class O, class C, class I> inline
00468 void dhorner( O& p, O& dp, C const * const c, unsigned n, const I& t, int s = 1 )
00469 {
00470 p = c[n-1], dp = 0.0;
00471 for ( int j = (n-2)*s; j >=0; dp=dp*t+p, p=p*t+c[j], j -= s ) ;
00472 };
00473
00474
00475 template<class A>
00476 void inverses( A * a, A * ea ) { for ( ;a != ea; *a++ = 1/(*a) ) ; };
00477
00478
00479 template<typename real_t>
00480 real_t min( real_t const * src, int sz, int st = 1 )
00481 {
00482 real_t r = src[0];
00483 for ( int p = st; p != sz*st; p += st ) if ( r > src[p] ) r = src[p];
00484 return r;
00485 };
00486
00487 template<typename real_t>
00488 real_t max( real_t const * src, int sz, int st = 1 )
00489 {
00490 real_t r = src[0];
00491 for ( int p = st; p != sz*st; p += st ) if ( r < src[p] ) r = src[p];
00492 return r;
00493 };
00494
00495
00496
00497
00498 };
00499
00500 }
00501
00502 #endif
