00001
00002
00003
00004
00005
00006 #ifndef realroot_SOLVE_SBDSLV_SUPPORT_HPP
00007 #define realroot_SOLVE_SBDSLV_SUPPORT_HPP
00008
00009 namespace mmx {
00010
00011 namespace realroot
00012 {
00013 template<class V>
00014 struct system_ctrl
00015 {
00016 typedef int sz_t;
00017 V& m_v;
00018 system_ctrl( V& v ): m_v(v) {};
00019 template< class _interval_ > static inline
00020 bool check( _interval_ const * const, int ) { return true; };
00021 template< class _interval_ > inline
00022 void output( _interval_ const * const dmns, sz_t nvars )
00023 {
00024
00025
00026
00027 for ( sz_t i = 0; i < nvars; i ++ )
00028 m_v.push_back( dmns[i].lower() ),
00029 m_v.push_back( dmns[i].upper() );
00030 };
00031 };
00032
00033 template< class real_t >
00034 unsigned clean_result( real_t * isols, int nvars, int nsols, const real_t& prec )
00035 {
00036 int a,c,n,i;
00037 real_t cd;
00038 if ( nsols == 0 ) return 0;
00039 for ( a = 0; a < nsols*nvars; isols[a] = (isols[2*a]+isols[2*a+1])/2.0, a ++ ) ;
00040
00041 for ( n = 1, c = 0, a = 0; a < nsols*nvars; a += nvars )
00042 {
00043 for ( cd = 0, i = 0; i < nvars; cd = std::max(cd,std::abs((isols+c)[i]-(isols+a)[i])), i ++ ) ;
00044 if ( cd < prec )
00045 {
00046 double fa = ((double)n)/(n+1);
00047 double fb = ((double)1.0)/(n+1);
00048 for ( int i = 0; i < nvars; i ++ ) (isols+c)[i] = fa*(isols+c)[i]+fb*(isols+a)[i];
00049 n ++;
00050 }
00051 else
00052 {
00053 c += nvars;
00054 std::copy(isols+a,isols+a+nvars,isols+c);
00055 n = 1;
00056 };
00057 };
00058 return (c+nvars)/nvars;
00059 };
00060 };
00061
00062 }
00063
00064 #endif //