|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
#include <solver_uv_interval_newton.hpp>
Definition at line 64 of file solver_uv_interval_newton.hpp.
| typedef INT::boundary_type BT |
Definition at line 67 of file solver_uv_interval_newton.hpp.
| typedef ring<CT, MonomialTensor>::Polynomial polynomial |
Definition at line 66 of file solver_uv_interval_newton.hpp.
| static INT bisection_iteration | ( | const polynomial & | p, |
| const INT | approx, | ||
| int & | status_ | ||
| ) | [inline, static] |
Definition at line 73 of file solver_uv_interval_newton.hpp.
References approx(), mmx::lower(), interval_newton< INT, CT >::median(), mmx::sign(), and mmx::upper().
Referenced by interval_newton< INT, CT >::newton_iteration().
{
BT x = median(approx);
INT lval, rval, mval;
lval = p(INT(mmx::lower(approx)));
rval = p(INT(mmx::upper(approx)));
mval = p(INT(x));
// std::cerr << lval << ' ' << rval << ' ' << mval << ' ' << x << std::endl;
// assert( mmx::sign(lval)*mmx::sign(rval) <= 0 );
status_ = 0;
if (mval == 0) return INT(x);
if (lval == 0) return INT(mmx::lower(approx));
if (rval == 0) return INT(mmx::upper(approx));
if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
status_ = -3; return approx;
}
if (mmx::sign(lval) == mmx::sign(rval)) {
status_ = -1; return approx;
}
if (mmx::sign(mval) == 0) return approx;
if (mmx::sign(lval) == mmx::sign(mval)) return INT(x,mmx::upper(approx));
return INT(mmx::lower(approx), x);
}
| static BT median | ( | const INT & | x | ) | [inline, static] |
Definition at line 70 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().
{ 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 112 of file solver_uv_interval_newton.hpp.
References approx(), assert, interval_newton< INT, CT >::bisection_iteration(), mmx::intersect(), mmx::lower(), interval_newton< INT, CT >::median(), mmx::sign(), mmx::upper(), and interval_newton< INT, CT >::width().
{
INT slope, testapprox, newapprox, ver;
BT oldw = width(approx);
BT neww;
int z;
#ifdef INEWTONDEBUG
int i = 0;
#endif
INT lval, rval;
lval = p(INT(mmx::lower(approx)));
rval = p(INT(mmx::upper(approx)));
#ifdef INEWTONDEBUG
std::cerr << "approx = " << approx << " lval = " << lval << " rval = " << rval << std::endl;
#endif
if (lval == 0) { status_ = 0; return INT(mmx::lower(approx)); }
if (rval == 0) { status_ = 0; return INT(mmx::upper(approx)); }
if (mmx::sign(lval) == 0 || mmx::sign(rval) == 0) {
status_ = -3; return approx;
}
status_ = 0;
slope = dp(approx);
#ifdef INEWTONDEBUG
std::cerr << "slope = " << slope << width(slope) << std::endl;
#endif
z = mmx::sign(slope);
if (mmx::sign(lval) == mmx::sign(rval)) {
if (z == 0) status_ = -1; else status_ = -2;
return approx;
}
#ifdef INEWTONDEBUG
std::cerr << "solving " << p << std::endl;
#endif
/* INTERVAL NEWTON MAIN LOOP */
while (1) {
// should we always do simple bisection with intervals of length >= 1 ?
// if (oldw >= 1 || z == 0) {
// slope contains zero
if (z == 0) {
approx = bisection_iteration(p, approx, status_);
neww = width(approx);
#ifdef INEWTONDEBUG
i++;
std::cerr << '(' << i << ") bisect: " << approx << ' ' << neww << std::endl;
#endif
if (status_ < 0) return approx;
if (neww < delta_) break;
assert(neww <= oldw);
if (neww == oldw) break;
oldw = neww;
} else {
// slope doesn't contain zero
BT x = median(approx);
if (x == mmx::lower(approx) || x == mmx::upper(approx)) break;
#ifdef INEWTONDEBUG
std::cerr << "mid = " << x << " mid_INT = " << INT(x) << std::endl;
std::cerr << "p(x) = " << p(INT(x)) << std::endl;
std::cerr << "p'(approx) = " << dp(approx) << std::endl;
std::cerr << "p/p' = " << INT(p(INT(x))/dp(approx)) << std::endl;
#endif
testapprox = INT(x) - INT(p(INT(x))/dp(approx)); // Newton approximation
#ifdef INEWTONDEBUG
i++;
#endif
if (mmx::intersect(newapprox, testapprox, approx)) {
#ifdef INEWTONDEBUG
std::cerr <<"new =" << newapprox << " test = " << testapprox << " old = " << approx << std::endl;
#endif
neww = width(newapprox);
// new approximation SHOULD contain the root
// if you want to verify, use the debugging code below
#ifdef INEWTONDEBUG
INT ver = p(newapprox);
if (ver != INT(0)) {
assert (mmx::sign(ver) == 0);
}
std::cerr << '(' << i << ") Newton: " << newapprox << ' ' << neww << std::endl;
#endif
//if (Float(2)*neww >= oldw) break;
if (neww >= oldw) break;
approx = newapprox;
if (neww < delta_) break;
oldw = neww;
} else {
// no convergence
status_ = -4;
return approx;
}
}
slope = dp(approx);
#ifdef INEWTONDEBUG
std::cerr << "slope = " << slope << width(slope) << std::endl;
#endif
z = mmx::sign(slope);
}
if (mmx::sign(dp(approx)) != 0) status_ = 1;
#ifdef INEWTONDEBUG
std::cerr << "iters = " << i << std::endl;
#endif
return approx;
}
| static bool singleton | ( | const INT & | x | ) | [inline, static] |
Definition at line 71 of file solver_uv_interval_newton.hpp.
References mmx::lower(), and mmx::upper().
{ return mmx::lower(x) == mmx::upper(x); }
| static BT width | ( | const INT & | x | ) | [inline, static] |
Definition at line 69 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().
{ return mmx::upper(x) - mmx::lower(x); }
