interval_newton< INT, CT > Struct Template Reference

#include <solver_uv_interval_newton.hpp>

List of all members.

Public Types

Static Public Member Functions

Detailed Description

template<class INT, class CT>
struct mmx::interval_newton< INT, CT >

Definition at line 65 of file solver_uv_interval_newton.hpp.

Member Typedef Documentation

typedef INT::boundary_type BT

Definition at line 68 of file solver_uv_interval_newton.hpp.

typedef ring<CT, MonomialTensor>::Polynomial polynomial

Definition at line 67 of file solver_uv_interval_newton.hpp.

Member Function Documentation

static INT bisection_iteration ( const polynomial &  p,
const INT  approx,
int &  status_ 
) [inline, static]

Definition at line 74 of file solver_uv_interval_newton.hpp.

References mmx::lower(), interval_newton< INT, CT >::median(), mmx::sign(), and mmx::upper().

Referenced by interval_newton< INT, CT >::newton_iteration().

00075       {
00076           BT x = median(approx);
00077           INT lval, rval, mval;
00078           lval = p(INT(mmx::lower(approx)));
00079           rval = p(INT(mmx::upper(approx)));
00080           mval = p(INT(x));
00081           //    std::cerr << lval << ' ' << rval << ' ' << mval << ' ' << x << std::endl;
00082           //    assert( mmx::sign(lval)*mmx::sign(rval) <= 0 );
00083                    
00084           status_ = 0;
00085           if (mval == 0) return INT(x);
00086           if (lval == 0) return INT(mmx::lower(approx)); 
00087           if (rval == 0) return INT(mmx::upper(approx));
00088                    
00089           if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
00090               status_ = -3; return approx;
00091           }
00092 
00093           if (mmx::sign(lval) == mmx::sign(rval)) {
00094               status_ = -1; return approx;
00095           }
00096           
00097           if (mmx::sign(mval) == 0) return approx;
00098           
00099           if (mmx::sign(lval) == mmx::sign(mval)) return INT(x,mmx::upper(approx));
00100           return INT(mmx::lower(approx), x);
00101       }

static BT median ( const INT &  x  )  [inline, static]

Definition at line 71 of file solver_uv_interval_newton.hpp.

References mmx::lower(), and interval_newton< INT, CT >::width().

Referenced by interval_newton< INT, CT >::bisection_iteration(), and interval_newton< INT, CT >::newton_iteration().

00071 { return mmx::lower(x) + width(x)/BT(2); }

static INT newton_iteration ( const polynomial &  p,
const polynomial &  dp,
INT  approx,
BT  delta_,
int &  status_ 
) [inline, static]

Definition at line 113 of file solver_uv_interval_newton.hpp.

References assert, interval_newton< INT, CT >::bisection_iteration(), mmx::intersect(), mmx::lower(), interval_newton< INT, CT >::median(), mmx::sign(), slope(), mmx::upper(), and interval_newton< INT, CT >::width().

00115       { 
00116           INT slope, testapprox, newapprox, ver;
00117           BT oldw = width(approx);
00118           BT neww;
00119           int z;
00120           
00121 #ifdef INEWTONDEBUG
00122           int i = 0;
00123 #endif
00124           
00125           INT lval, rval;
00126           
00127           lval = p(INT(mmx::lower(approx)));
00128           rval = p(INT(mmx::upper(approx)));
00129 #ifdef INEWTONDEBUG
00130           std::cerr << "approx = " << approx << " lval = " << lval << " rval = " << rval << std::endl;
00131 #endif
00132           if (lval == 0) { status_ = 0; return INT(mmx::lower(approx)); }
00133           if (rval == 0) { status_ = 0; return INT(mmx::upper(approx)); }
00134           if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
00135               status_ = -3; return approx;
00136           }
00137           
00138           status_ = 0;
00139           slope = dp(approx);
00140           
00141 #ifdef INEWTONDEBUG
00142           std::cerr << "slope = " << slope << width(slope) << std::endl;
00143 #endif
00144 
00145           z = mmx::sign(slope);
00146           if (mmx::sign(lval) == mmx::sign(rval)) {
00147               if (z == 0) status_ = -1; else status_ = -2;
00148               return approx;
00149           }
00150 
00151 #ifdef INEWTONDEBUG
00152           std::cerr << "solving " << p << std::endl;
00153 #endif
00154 
00155           /* INTERVAL NEWTON MAIN LOOP */
00156           while (1) {
00157               // should we always do simple bisection with intervals of length >= 1 ?
00158               // if (oldw >= 1 || z == 0) {
00159               
00160               // slope contains zero
00161               if (z == 0) {
00162                   approx = bisection_iteration(p, approx, status_);
00163                   neww = width(approx);
00164 #ifdef INEWTONDEBUG
00165                   i++;
00166                   std::cerr << '(' << i << ") bisect: " << approx << ' ' << neww << std::endl;
00167 #endif
00168                   if (status_ < 0) return approx;
00169                   if (neww < delta_) break;
00170                   assert(neww <= oldw);
00171                   if (neww == oldw) break;
00172                   oldw = neww;
00173                   
00174               } else {
00175               // slope doesn't contain zero
00176 
00177                   BT x = median(approx);
00178                   if (x == mmx::lower(approx) || x == mmx::upper(approx)) break;
00179 #ifdef INEWTONDEBUG
00180                   std::cerr << "mid = " << x << " mid_INT = " << INT(x) << std::endl;
00181                   std::cerr << "p(x) = " << p(INT(x)) << std::endl;
00182                   std::cerr << "p'(approx) = " << dp(approx) << std::endl;
00183                   std::cerr << "p/p' = " << INT(p(INT(x))/dp(approx)) << std::endl;
00184 #endif
00185                   testapprox = INT(x) - INT(p(INT(x))/dp(approx)); // Newton approximation
00186 
00187 #ifdef INEWTONDEBUG
00188                   i++;
00189 #endif
00190                   if (mmx::intersect(newapprox, testapprox, approx)) {
00191 #ifdef INEWTONDEBUG
00192                       std::cerr <<"new =" << newapprox << " test = " << testapprox << " old = " << approx << std::endl;
00193 #endif
00194                       neww = width(newapprox);
00195                       
00196                       // new approximation SHOULD contain the root
00197                       // if you want to verify, use the debugging code below
00198 #ifdef INEWTONDEBUG
00199                       INT ver = p(newapprox);
00200                       if (ver != INT(0)) {
00201                           assert (mmx::sign(ver) == 0);
00202                       }
00203                       
00204                       std::cerr << '(' << i << ") Newton: " << newapprox << ' ' << neww << std::endl;
00205 #endif
00206                       //if (Float(2)*neww >= oldw) break;
00207                       if (neww >= oldw) break;
00208                       approx = newapprox;
00209                       if (neww < delta_) break;
00210                       oldw = neww;
00211                   } else {
00212                       // no convergence
00213                       status_ = -4;
00214                       return approx;
00215                   }
00216               }
00217               slope = dp(approx);
00218 #ifdef INEWTONDEBUG
00219               std::cerr << "slope = " << slope << width(slope) << std::endl;
00220 #endif
00221               z = mmx::sign(slope);
00222           }
00223           if (mmx::sign(dp(approx)) != 0) status_ = 1;
00224 #ifdef INEWTONDEBUG
00225           std::cerr << "iters = " << i << std::endl;
00226 #endif
00227           return approx;
00228       }

static bool singleton ( const INT &  x  )  [inline, static]

Definition at line 72 of file solver_uv_interval_newton.hpp.

References mmx::lower(), and mmx::upper().

00072 { return mmx::lower(x) == mmx::upper(x); }

static BT width ( const INT &  x  )  [inline, static]

Definition at line 70 of file solver_uv_interval_newton.hpp.

References mmx::lower(), and mmx::upper().

Referenced by interval_newton< INT, CT >::median(), and interval_newton< INT, CT >::newton_iteration().

00070 { return mmx::upper(x) - mmx::lower(x); }

The documentation for this struct was generated from the following file:
Generated on 6 Dec 2012 for realroot by  doxygen 1.6.1