- Merge with the static filters.

- Added CGAL_assertion(Interval_nt_advanced::want_exceptions) for the dynamic.
This commit is contained in:
Sylvain Pion 1999-08-05 18:29:01 +00:00
parent ad1bb0b0bd
commit 8f77b04a95
8 changed files with 6150 additions and 149 deletions

View File

@ -26,6 +26,192 @@
CGAL_BEGIN_NAMESPACE
inline
Oriented_side
power_testC2_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pwt,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qwt,
const Static_filter_error &rx,
const Static_filter_error &ry,
const Static_filter_error &rwt,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &twt,
double & epsilon_0)
{
typedef Static_filter_error FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = square(dpx) + square(dpy) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = square(dqx) + square(dqy) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = square(drx) + square(dry) - rwt + twt;
return Oriented_side(sign_of_determinant3x3_SAF(dpx, dpy, dpz,
dqx, dqy, dqz,
drx, dry, drz,
epsilon_0));
}
inline
Oriented_side
power_testC2_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pwt,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qwt,
const Restricted_double &rx,
const Restricted_double &ry,
const Restricted_double &rwt,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &twt,
const double & epsilon_0)
{
typedef Restricted_double FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = square(dpx) + square(dpy) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = square(dqx) + square(dqy) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = square(drx) + square(dry) - rwt + twt;
return Oriented_side(sign_of_determinant3x3_SAF(dpx, dpy, dpz,
dqx, dqy, dqz,
drx, dry, drz,
epsilon_0));
}
inline
Oriented_side
power_testC2(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &rx,
const Static_adaptatif_filter &ry,
const Static_adaptatif_filter &rwt,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(rx.value()) > SAF_bound ||
fabs(ry.value()) > SAF_bound ||
fabs(rwt.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(rx.value()));
SAF_bound = std::max(SAF_bound, fabs(ry.value()));
SAF_bound = std::max(SAF_bound, fabs(rwt.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testC2_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testC2_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pwt.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qwt.value()),
Restricted_double(rx.value()),
Restricted_double(ry.value()),
Restricted_double(rwt.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(twt.value()),
SAF_epsilon_0);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testC2(
px.exact(),
py.exact(),
pwt.exact(),
qx.exact(),
qy.exact(),
qwt.exact(),
rx.exact(),
ry.exact(),
rwt.exact(),
tx.exact(),
ty.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -45,6 +231,7 @@ power_testC2(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &ty,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{
@ -83,6 +270,189 @@ power_testC2(
}
}
inline
Oriented_side
power_testC2_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pwt,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qwt,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &twt,
double & epsilon_0,
double & epsilon_1,
double & epsilon_2,
double & epsilon_3)
{
typedef Static_filter_error FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = square(dpx) + square(dpy) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = square(dqx) + square(dqy) - qwt + twt;
Comparison_result cmpx = CGAL::compare_SAF(px, qx,
epsilon_0);
if (cmpx != EQUAL)
return Oriented_side(cmpx * sign_of_determinant2x2_SAF(dpx, dpz, dqx, dqz,
epsilon_1));
Comparison_result cmpy = CGAL::compare_SAF(py, qy,
epsilon_2);
return Oriented_side(cmpy * sign_of_determinant2x2_SAF(dpy, dpz, dqy, dqz,
epsilon_3));
}
inline
Oriented_side
power_testC2_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pwt,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qwt,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &twt,
const double & epsilon_0,
const double & epsilon_1,
const double & epsilon_2,
const double & epsilon_3)
{
typedef Restricted_double FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = square(dpx) + square(dpy) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = square(dqx) + square(dqy) - qwt + twt;
Comparison_result cmpx = CGAL::compare_SAF(px, qx,
epsilon_0);
if (cmpx != EQUAL)
return Oriented_side(cmpx * sign_of_determinant2x2_SAF(dpx, dpz, dqx, dqz,
epsilon_1));
Comparison_result cmpy = CGAL::compare_SAF(py, qy,
epsilon_2);
return Oriented_side(cmpy * sign_of_determinant2x2_SAF(dpy, dpz, dqy, dqz,
epsilon_3));
}
inline
Oriented_side
power_testC2(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
static double SAF_epsilon_1;
static double SAF_epsilon_2;
static double SAF_epsilon_3;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testC2_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testC2_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pwt.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qwt.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(twt.value()),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testC2(
px.exact(),
py.exact(),
pwt.exact(),
qx.exact(),
qy.exact(),
qwt.exact(),
tx.exact(),
ty.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -99,6 +469,7 @@ power_testC2(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &ty,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{

View File

@ -26,6 +26,268 @@
CGAL_BEGIN_NAMESPACE
inline
Oriented_side
power_testC3_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pz,
const Static_filter_error &pwt,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qz,
const Static_filter_error &qwt,
const Static_filter_error &rx,
const Static_filter_error &ry,
const Static_filter_error &rz,
const Static_filter_error &rwt,
const Static_filter_error &sx,
const Static_filter_error &sy,
const Static_filter_error &sz,
const Static_filter_error &swt,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &tz,
const Static_filter_error &twt,
double & epsilon_0)
{
typedef Static_filter_error FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = rz - tz;
FT drt = square(drx) + square(dry) + square(drz) - rwt + twt;
FT dsx = sx - tx;
FT dsy = sy - ty;
FT dsz = sz - tz;
FT dst = square(dsx) + square(dsy) + square(dsz) - swt + twt;
return Oriented_side( - sign_of_determinant4x4_SAF(dpx, dpy, dpz, dpt,
dqx, dqy, dqz, dqt,
drx, dry, drz, drt,
dsx, dsy, dsz, dst,
epsilon_0));
}
inline
Oriented_side
power_testC3_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pz,
const Restricted_double &pwt,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qz,
const Restricted_double &qwt,
const Restricted_double &rx,
const Restricted_double &ry,
const Restricted_double &rz,
const Restricted_double &rwt,
const Restricted_double &sx,
const Restricted_double &sy,
const Restricted_double &sz,
const Restricted_double &swt,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &tz,
const Restricted_double &twt,
const double & epsilon_0)
{
typedef Restricted_double FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = rz - tz;
FT drt = square(drx) + square(dry) + square(drz) - rwt + twt;
FT dsx = sx - tx;
FT dsy = sy - ty;
FT dsz = sz - tz;
FT dst = square(dsx) + square(dsy) + square(dsz) - swt + twt;
return Oriented_side( - sign_of_determinant4x4_SAF(dpx, dpy, dpz, dpt,
dqx, dqy, dqz, dqt,
drx, dry, drz, drt,
dsx, dsy, dsz, dst,
epsilon_0));
}
inline
Oriented_side
power_testC3(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pz,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qz,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &rx,
const Static_adaptatif_filter &ry,
const Static_adaptatif_filter &rz,
const Static_adaptatif_filter &rwt,
const Static_adaptatif_filter &sx,
const Static_adaptatif_filter &sy,
const Static_adaptatif_filter &sz,
const Static_adaptatif_filter &swt,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &tz,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pz.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qz.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(rx.value()) > SAF_bound ||
fabs(ry.value()) > SAF_bound ||
fabs(rz.value()) > SAF_bound ||
fabs(rwt.value()) > SAF_bound ||
fabs(sx.value()) > SAF_bound ||
fabs(sy.value()) > SAF_bound ||
fabs(sz.value()) > SAF_bound ||
fabs(swt.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(tz.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pz.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qz.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(rx.value()));
SAF_bound = std::max(SAF_bound, fabs(ry.value()));
SAF_bound = std::max(SAF_bound, fabs(rz.value()));
SAF_bound = std::max(SAF_bound, fabs(rwt.value()));
SAF_bound = std::max(SAF_bound, fabs(sx.value()));
SAF_bound = std::max(SAF_bound, fabs(sy.value()));
SAF_bound = std::max(SAF_bound, fabs(sz.value()));
SAF_bound = std::max(SAF_bound, fabs(swt.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(tz.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testC3_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testC3_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pz.value()),
Restricted_double(pwt.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qz.value()),
Restricted_double(qwt.value()),
Restricted_double(rx.value()),
Restricted_double(ry.value()),
Restricted_double(rz.value()),
Restricted_double(rwt.value()),
Restricted_double(sx.value()),
Restricted_double(sy.value()),
Restricted_double(sz.value()),
Restricted_double(swt.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(tz.value()),
Restricted_double(twt.value()),
SAF_epsilon_0);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testC3(
px.exact(),
py.exact(),
pz.exact(),
pwt.exact(),
qx.exact(),
qy.exact(),
qz.exact(),
qwt.exact(),
rx.exact(),
ry.exact(),
rz.exact(),
rwt.exact(),
sx.exact(),
sy.exact(),
sz.exact(),
swt.exact(),
tx.exact(),
ty.exact(),
tz.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -108,6 +370,301 @@ power_testC3(
}
}
inline
Oriented_side
power_testC3_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pz,
const Static_filter_error &pwt,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qz,
const Static_filter_error &qwt,
const Static_filter_error &rx,
const Static_filter_error &ry,
const Static_filter_error &rz,
const Static_filter_error &rwt,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &tz,
const Static_filter_error &twt,
double & epsilon_0,
double & epsilon_1,
double & epsilon_2,
double & epsilon_3,
double & epsilon_4,
double & epsilon_5)
{
typedef Static_filter_error FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = rz - tz;
FT drt = square(drx) + square(dry) + square(drz) - rwt + twt;
Sign cmp;
cmp = sign_of_determinant3x3_SAF(dpx, dpy, dpt,
dqx, dqy, dqt,
drx, dry, drt,
epsilon_0);
if (cmp != ZERO)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, py-ry,
qx-rx, qy-ry,
epsilon_1));
cmp = sign_of_determinant3x3_SAF(dpx, dpz, dpt,
dqx, dqz, dqt,
drx, drz, drt,
epsilon_2);
if (cmp != ZERO)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, pz-rz,
qx-rx, qz-rz,
epsilon_3));
cmp = sign_of_determinant3x3_SAF(dpy, dpz, dpt,
dqy, dqz, dqt,
dry, drz, drt,
epsilon_4);
return Oriented_side(cmp * sign_of_determinant2x2_SAF(py-ry, pz-rz,
qy-ry, qz-rz,
epsilon_5));
}
inline
Oriented_side
power_testC3_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pz,
const Restricted_double &pwt,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qz,
const Restricted_double &qwt,
const Restricted_double &rx,
const Restricted_double &ry,
const Restricted_double &rz,
const Restricted_double &rwt,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &tz,
const Restricted_double &twt,
const double & epsilon_0,
const double & epsilon_1,
const double & epsilon_2,
const double & epsilon_3,
const double & epsilon_4,
const double & epsilon_5)
{
typedef Restricted_double FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
FT drx = rx - tx;
FT dry = ry - ty;
FT drz = rz - tz;
FT drt = square(drx) + square(dry) + square(drz) - rwt + twt;
Sign cmp;
cmp = sign_of_determinant3x3_SAF(dpx, dpy, dpt,
dqx, dqy, dqt,
drx, dry, drt,
epsilon_0);
if (cmp != ZERO)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, py-ry,
qx-rx, qy-ry,
epsilon_1));
cmp = sign_of_determinant3x3_SAF(dpx, dpz, dpt,
dqx, dqz, dqt,
drx, drz, drt,
epsilon_2);
if (cmp != ZERO)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, pz-rz,
qx-rx, qz-rz,
epsilon_3));
cmp = sign_of_determinant3x3_SAF(dpy, dpz, dpt,
dqy, dqz, dqt,
dry, drz, drt,
epsilon_4);
return Oriented_side(cmp * sign_of_determinant2x2_SAF(py-ry, pz-rz,
qy-ry, qz-rz,
epsilon_5));
}
inline
Oriented_side
power_testC3(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pz,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qz,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &rx,
const Static_adaptatif_filter &ry,
const Static_adaptatif_filter &rz,
const Static_adaptatif_filter &rwt,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &tz,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
static double SAF_epsilon_1;
static double SAF_epsilon_2;
static double SAF_epsilon_3;
static double SAF_epsilon_4;
static double SAF_epsilon_5;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pz.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qz.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(rx.value()) > SAF_bound ||
fabs(ry.value()) > SAF_bound ||
fabs(rz.value()) > SAF_bound ||
fabs(rwt.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(tz.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pz.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qz.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(rx.value()));
SAF_bound = std::max(SAF_bound, fabs(ry.value()));
SAF_bound = std::max(SAF_bound, fabs(rz.value()));
SAF_bound = std::max(SAF_bound, fabs(rwt.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(tz.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testC3_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3,
SAF_epsilon_4,
SAF_epsilon_5);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testC3_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pz.value()),
Restricted_double(pwt.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qz.value()),
Restricted_double(qwt.value()),
Restricted_double(rx.value()),
Restricted_double(ry.value()),
Restricted_double(rz.value()),
Restricted_double(rwt.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(tz.value()),
Restricted_double(twt.value()),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3,
SAF_epsilon_4,
SAF_epsilon_5);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testC3(
px.exact(),
py.exact(),
pz.exact(),
pwt.exact(),
qx.exact(),
qy.exact(),
qz.exact(),
qwt.exact(),
rx.exact(),
ry.exact(),
rz.exact(),
rwt.exact(),
tx.exact(),
ty.exact(),
tz.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -178,6 +735,243 @@ power_testC3(
}
}
inline
Oriented_side
power_testC3_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pz,
const Static_filter_error &pwt,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qz,
const Static_filter_error &qwt,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &tz,
const Static_filter_error &twt,
double & epsilon_0,
double & epsilon_1,
double & epsilon_2,
double & epsilon_3,
double & epsilon_4,
double & epsilon_5)
{
typedef Static_filter_error FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
Comparison_result cmp;
cmp = CGAL::compare_SAF(px, qx,
epsilon_0);
if (cmp != EQUAL)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpx, dpt, dqx, dqt,
epsilon_1));
cmp = CGAL::compare_SAF(py, qy,
epsilon_2);
if (cmp != EQUAL)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpy, dpt, dqy, dqt,
epsilon_3));
cmp = CGAL::compare_SAF(pz, qz,
epsilon_4);
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpz, dpt, dqz, dqt,
epsilon_5));
}
inline
Oriented_side
power_testC3_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pz,
const Restricted_double &pwt,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qz,
const Restricted_double &qwt,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &tz,
const Restricted_double &twt,
const double & epsilon_0,
const double & epsilon_1,
const double & epsilon_2,
const double & epsilon_3,
const double & epsilon_4,
const double & epsilon_5)
{
typedef Restricted_double FT;
FT dpx = px - tx;
FT dpy = py - ty;
FT dpz = pz - tz;
FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt;
FT dqx = qx - tx;
FT dqy = qy - ty;
FT dqz = qz - tz;
FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt;
Comparison_result cmp;
cmp = CGAL::compare_SAF(px, qx,
epsilon_0);
if (cmp != EQUAL)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpx, dpt, dqx, dqt,
epsilon_1));
cmp = CGAL::compare_SAF(py, qy,
epsilon_2);
if (cmp != EQUAL)
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpy, dpt, dqy, dqt,
epsilon_3));
cmp = CGAL::compare_SAF(pz, qz,
epsilon_4);
return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpz, dpt, dqz, dqt,
epsilon_5));
}
inline
Oriented_side
power_testC3(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pz,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qz,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &tz,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
static double SAF_epsilon_1;
static double SAF_epsilon_2;
static double SAF_epsilon_3;
static double SAF_epsilon_4;
static double SAF_epsilon_5;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pz.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qz.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(tz.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pz.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qz.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(tz.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testC3_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3,
SAF_epsilon_4,
SAF_epsilon_5);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testC3_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pz.value()),
Restricted_double(pwt.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qz.value()),
Restricted_double(qwt.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(tz.value()),
Restricted_double(twt.value()),
SAF_epsilon_0,
SAF_epsilon_1,
SAF_epsilon_2,
SAF_epsilon_3,
SAF_epsilon_4,
SAF_epsilon_5);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testC3(
px.exact(),
py.exact(),
pz.exact(),
pwt.exact(),
qx.exact(),
qy.exact(),
qz.exact(),
qwt.exact(),
tx.exact(),
ty.exact(),
tz.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif

View File

@ -26,6 +26,240 @@
CGAL_BEGIN_NAMESPACE
inline
Oriented_side
power_testH2_SAF(
const Static_filter_error &phx,
const Static_filter_error &phy,
const Static_filter_error &phw,
const Static_filter_error &pwt,
const Static_filter_error &qhx,
const Static_filter_error &qhy,
const Static_filter_error &qhw,
const Static_filter_error &qwt,
const Static_filter_error &rhx,
const Static_filter_error &rhy,
const Static_filter_error &rhw,
const Static_filter_error &rwt,
const Static_filter_error &thx,
const Static_filter_error &thy,
const Static_filter_error &thw,
const Static_filter_error &twt,
double & epsilon_0)
{
typedef Static_filter_error RT;
RT dphx = phx*phw;
RT dphy = phy*phw;
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) - pwt*dphw;
RT dqhx = qhx*qhw;
RT dqhy = qhy*qhw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) - qwt*dqhw;
RT drhx = rhx*rhw;
RT drhy = rhy*rhw;
RT drhw = square(rhw);
RT drz = square(rhx) + square(rhy) - rwt*drhw;
RT dthx = thx*thw;
RT dthy = thy*thw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) - twt*dthw;
return Oriented_side(sign_of_determinant4x4_SAF(dphx, dphy, dpz, dphw,
dqhx, dqhy, dqz, dqhw,
drhx, drhy, drz, drhw,
dthx, dthy, dtz, dthw,
epsilon_0));
}
inline
Oriented_side
power_testH2_SAF(
const Restricted_double &phx,
const Restricted_double &phy,
const Restricted_double &phw,
const Restricted_double &pwt,
const Restricted_double &qhx,
const Restricted_double &qhy,
const Restricted_double &qhw,
const Restricted_double &qwt,
const Restricted_double &rhx,
const Restricted_double &rhy,
const Restricted_double &rhw,
const Restricted_double &rwt,
const Restricted_double &thx,
const Restricted_double &thy,
const Restricted_double &thw,
const Restricted_double &twt,
const double & epsilon_0)
{
typedef Restricted_double RT;
RT dphx = phx*phw;
RT dphy = phy*phw;
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) - pwt*dphw;
RT dqhx = qhx*qhw;
RT dqhy = qhy*qhw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) - qwt*dqhw;
RT drhx = rhx*rhw;
RT drhy = rhy*rhw;
RT drhw = square(rhw);
RT drz = square(rhx) + square(rhy) - rwt*drhw;
RT dthx = thx*thw;
RT dthy = thy*thw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) - twt*dthw;
return Oriented_side(sign_of_determinant4x4_SAF(dphx, dphy, dpz, dphw,
dqhx, dqhy, dqz, dqhw,
drhx, drhy, drz, drhw,
dthx, dthy, dtz, dthw,
epsilon_0));
}
inline
Oriented_side
power_testH2(
const Static_adaptatif_filter &phx,
const Static_adaptatif_filter &phy,
const Static_adaptatif_filter &phw,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qhx,
const Static_adaptatif_filter &qhy,
const Static_adaptatif_filter &qhw,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &rhx,
const Static_adaptatif_filter &rhy,
const Static_adaptatif_filter &rhw,
const Static_adaptatif_filter &rwt,
const Static_adaptatif_filter &thx,
const Static_adaptatif_filter &thy,
const Static_adaptatif_filter &thw,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(phx.value()) > SAF_bound ||
fabs(phy.value()) > SAF_bound ||
fabs(phw.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qhx.value()) > SAF_bound ||
fabs(qhy.value()) > SAF_bound ||
fabs(qhw.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(rhx.value()) > SAF_bound ||
fabs(rhy.value()) > SAF_bound ||
fabs(rhw.value()) > SAF_bound ||
fabs(rwt.value()) > SAF_bound ||
fabs(thx.value()) > SAF_bound ||
fabs(thy.value()) > SAF_bound ||
fabs(thw.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(phx.value()));
SAF_bound = std::max(SAF_bound, fabs(phy.value()));
SAF_bound = std::max(SAF_bound, fabs(phw.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qhx.value()));
SAF_bound = std::max(SAF_bound, fabs(qhy.value()));
SAF_bound = std::max(SAF_bound, fabs(qhw.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(rhx.value()));
SAF_bound = std::max(SAF_bound, fabs(rhy.value()));
SAF_bound = std::max(SAF_bound, fabs(rhw.value()));
SAF_bound = std::max(SAF_bound, fabs(rwt.value()));
SAF_bound = std::max(SAF_bound, fabs(thx.value()));
SAF_bound = std::max(SAF_bound, fabs(thy.value()));
SAF_bound = std::max(SAF_bound, fabs(thw.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testH2_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testH2_SAF(
Restricted_double(phx.value()),
Restricted_double(phy.value()),
Restricted_double(phw.value()),
Restricted_double(pwt.value()),
Restricted_double(qhx.value()),
Restricted_double(qhy.value()),
Restricted_double(qhw.value()),
Restricted_double(qwt.value()),
Restricted_double(rhx.value()),
Restricted_double(rhy.value()),
Restricted_double(rhw.value()),
Restricted_double(rwt.value()),
Restricted_double(thx.value()),
Restricted_double(thy.value()),
Restricted_double(thw.value()),
Restricted_double(twt.value()),
SAF_epsilon_0);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testH2(
phx.exact(),
phy.exact(),
phw.exact(),
pwt.exact(),
qhx.exact(),
qhy.exact(),
qhw.exact(),
qwt.exact(),
rhx.exact(),
rhy.exact(),
rhw.exact(),
rwt.exact(),
thx.exact(),
thy.exact(),
thw.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -49,6 +283,7 @@ power_testH2(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{
@ -95,6 +330,227 @@ power_testH2(
}
}
inline
Oriented_side
power_testH2_SAF(
const Static_filter_error &phx,
const Static_filter_error &phy,
const Static_filter_error &phw,
const Static_filter_error &pwt,
const Static_filter_error &qhx,
const Static_filter_error &qhy,
const Static_filter_error &qhw,
const Static_filter_error &qwt,
const Static_filter_error &thx,
const Static_filter_error &thy,
const Static_filter_error &thw,
const Static_filter_error &twt,
double & epsilon_0,
double & epsilon_1)
{
typedef Static_filter_error RT;
RT pa, qa, ta;
if (phx * qhw != qhx * phw )
{
pa = phx*phw;
qa = qhx*qhw;
ta = thx*thw;
}
else
{
pa = phy*phw;
qa = qhy*qhw;
ta = thy*thw;
}
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) - pwt*dphw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) - qwt*dqhw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) - twt*dthw;
return Oriented_side(CGAL::compare_SAF(pa, qa,
epsilon_0) *
sign_of_determinant3x3_SAF(pa, dpz, dphw,
qa, dqz, dqhw,
ta, dtz, dthw,
epsilon_1));
}
inline
Oriented_side
power_testH2_SAF(
const Restricted_double &phx,
const Restricted_double &phy,
const Restricted_double &phw,
const Restricted_double &pwt,
const Restricted_double &qhx,
const Restricted_double &qhy,
const Restricted_double &qhw,
const Restricted_double &qwt,
const Restricted_double &thx,
const Restricted_double &thy,
const Restricted_double &thw,
const Restricted_double &twt,
const double & epsilon_0,
const double & epsilon_1)
{
typedef Restricted_double RT;
RT pa, qa, ta;
if (phx * qhw != qhx * phw )
{
pa = phx*phw;
qa = qhx*qhw;
ta = thx*thw;
}
else
{
pa = phy*phw;
qa = qhy*qhw;
ta = thy*thw;
}
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) - pwt*dphw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) - qwt*dqhw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) - twt*dthw;
return Oriented_side(CGAL::compare_SAF(pa, qa,
epsilon_0) *
sign_of_determinant3x3_SAF(pa, dpz, dphw,
qa, dqz, dqhw,
ta, dtz, dthw,
epsilon_1));
}
inline
Oriented_side
power_testH2(
const Static_adaptatif_filter &phx,
const Static_adaptatif_filter &phy,
const Static_adaptatif_filter &phw,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qhx,
const Static_adaptatif_filter &qhy,
const Static_adaptatif_filter &qhw,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &thx,
const Static_adaptatif_filter &thy,
const Static_adaptatif_filter &thw,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
static double SAF_epsilon_1;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(phx.value()) > SAF_bound ||
fabs(phy.value()) > SAF_bound ||
fabs(phw.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qhx.value()) > SAF_bound ||
fabs(qhy.value()) > SAF_bound ||
fabs(qhw.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(thx.value()) > SAF_bound ||
fabs(thy.value()) > SAF_bound ||
fabs(thw.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(phx.value()));
SAF_bound = std::max(SAF_bound, fabs(phy.value()));
SAF_bound = std::max(SAF_bound, fabs(phw.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qhx.value()));
SAF_bound = std::max(SAF_bound, fabs(qhy.value()));
SAF_bound = std::max(SAF_bound, fabs(qhw.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(thx.value()));
SAF_bound = std::max(SAF_bound, fabs(thy.value()));
SAF_bound = std::max(SAF_bound, fabs(thw.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testH2_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0,
SAF_epsilon_1);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testH2_SAF(
Restricted_double(phx.value()),
Restricted_double(phy.value()),
Restricted_double(phw.value()),
Restricted_double(pwt.value()),
Restricted_double(qhx.value()),
Restricted_double(qhy.value()),
Restricted_double(qhw.value()),
Restricted_double(qwt.value()),
Restricted_double(thx.value()),
Restricted_double(thy.value()),
Restricted_double(thw.value()),
Restricted_double(twt.value()),
SAF_epsilon_0,
SAF_epsilon_1);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testH2(
phx.exact(),
phy.exact(),
phw.exact(),
pwt.exact(),
qhx.exact(),
qhy.exact(),
qhw.exact(),
qwt.exact(),
thx.exact(),
thy.exact(),
thw.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -114,6 +570,7 @@ power_testH2(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{

View File

@ -26,6 +26,334 @@
CGAL_BEGIN_NAMESPACE
inline
Oriented_side
power_testH3_SAF(
const Static_filter_error &phx,
const Static_filter_error &phy,
const Static_filter_error &phz,
const Static_filter_error &phw,
const Static_filter_error &pwt,
const Static_filter_error &qhx,
const Static_filter_error &qhy,
const Static_filter_error &qhz,
const Static_filter_error &qhw,
const Static_filter_error &qwt,
const Static_filter_error &rhx,
const Static_filter_error &rhy,
const Static_filter_error &rhz,
const Static_filter_error &rhw,
const Static_filter_error &rwt,
const Static_filter_error &shx,
const Static_filter_error &shy,
const Static_filter_error &shz,
const Static_filter_error &shw,
const Static_filter_error &swt,
const Static_filter_error &thx,
const Static_filter_error &thy,
const Static_filter_error &thz,
const Static_filter_error &thw,
const Static_filter_error &twt,
double & epsilon_0)
{
typedef Static_filter_error RT;
RT dphx = phx*phw;
RT dphy = phy*phw;
RT dphz = phz*phw;
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) + square(phz) - pwt*dphw;
RT dqhx = qhx*qhw;
RT dqhy = qhy*qhw;
RT dqhz = qhz*qhw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) + square(qhz) - qwt*dqhw;
RT drhx = rhx*rhw;
RT drhy = rhy*rhw;
RT drhz = rhz*rhw;
RT drhw = square(rhw);
RT drz = square(rhx) + square(rhy) + square(rhz) - rwt*drhw;
RT dshx = shx*shw;
RT dshy = shy*shw;
RT dshz = shz*shw;
RT dshw = square(shw);
RT dsz = square(shx) + square(shy) + square(shz) - swt*dshw;
RT dthx = thx*thw;
RT dthy = thy*thw;
RT dthz = thz*thw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) + square(thz) - twt*dthw;
return Oriented_side(- sign_of_determinant5x5_SAF(dphx, dphy, dphz, dpz, dphw,
dqhx, dqhy, dqhz, dqz, dqhw,
drhx, drhy, drhz, drz, drhw,
dshx, dshy, dshz, dsz, dshw,
dthx, dthy, dthz, dtz, dthw,
epsilon_0));
}
inline
Oriented_side
power_testH3_SAF(
const Restricted_double &phx,
const Restricted_double &phy,
const Restricted_double &phz,
const Restricted_double &phw,
const Restricted_double &pwt,
const Restricted_double &qhx,
const Restricted_double &qhy,
const Restricted_double &qhz,
const Restricted_double &qhw,
const Restricted_double &qwt,
const Restricted_double &rhx,
const Restricted_double &rhy,
const Restricted_double &rhz,
const Restricted_double &rhw,
const Restricted_double &rwt,
const Restricted_double &shx,
const Restricted_double &shy,
const Restricted_double &shz,
const Restricted_double &shw,
const Restricted_double &swt,
const Restricted_double &thx,
const Restricted_double &thy,
const Restricted_double &thz,
const Restricted_double &thw,
const Restricted_double &twt,
const double & epsilon_0)
{
typedef Restricted_double RT;
RT dphx = phx*phw;
RT dphy = phy*phw;
RT dphz = phz*phw;
RT dphw = square(phw);
RT dpz = square(phx) + square(phy) + square(phz) - pwt*dphw;
RT dqhx = qhx*qhw;
RT dqhy = qhy*qhw;
RT dqhz = qhz*qhw;
RT dqhw = square(qhw);
RT dqz = square(qhx) + square(qhy) + square(qhz) - qwt*dqhw;
RT drhx = rhx*rhw;
RT drhy = rhy*rhw;
RT drhz = rhz*rhw;
RT drhw = square(rhw);
RT drz = square(rhx) + square(rhy) + square(rhz) - rwt*drhw;
RT dshx = shx*shw;
RT dshy = shy*shw;
RT dshz = shz*shw;
RT dshw = square(shw);
RT dsz = square(shx) + square(shy) + square(shz) - swt*dshw;
RT dthx = thx*thw;
RT dthy = thy*thw;
RT dthz = thz*thw;
RT dthw = square(thw);
RT dtz = square(thx) + square(thy) + square(thz) - twt*dthw;
return Oriented_side(- sign_of_determinant5x5_SAF(dphx, dphy, dphz, dpz, dphw,
dqhx, dqhy, dqhz, dqz, dqhw,
drhx, drhy, drhz, drz, drhw,
dshx, dshy, dshz, dsz, dshw,
dthx, dthy, dthz, dtz, dthw,
epsilon_0));
}
inline
Oriented_side
power_testH3(
const Static_adaptatif_filter &phx,
const Static_adaptatif_filter &phy,
const Static_adaptatif_filter &phz,
const Static_adaptatif_filter &phw,
const Static_adaptatif_filter &pwt,
const Static_adaptatif_filter &qhx,
const Static_adaptatif_filter &qhy,
const Static_adaptatif_filter &qhz,
const Static_adaptatif_filter &qhw,
const Static_adaptatif_filter &qwt,
const Static_adaptatif_filter &rhx,
const Static_adaptatif_filter &rhy,
const Static_adaptatif_filter &rhz,
const Static_adaptatif_filter &rhw,
const Static_adaptatif_filter &rwt,
const Static_adaptatif_filter &shx,
const Static_adaptatif_filter &shy,
const Static_adaptatif_filter &shz,
const Static_adaptatif_filter &shw,
const Static_adaptatif_filter &swt,
const Static_adaptatif_filter &thx,
const Static_adaptatif_filter &thy,
const Static_adaptatif_filter &thz,
const Static_adaptatif_filter &thw,
const Static_adaptatif_filter &twt)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(phx.value()) > SAF_bound ||
fabs(phy.value()) > SAF_bound ||
fabs(phz.value()) > SAF_bound ||
fabs(phw.value()) > SAF_bound ||
fabs(pwt.value()) > SAF_bound ||
fabs(qhx.value()) > SAF_bound ||
fabs(qhy.value()) > SAF_bound ||
fabs(qhz.value()) > SAF_bound ||
fabs(qhw.value()) > SAF_bound ||
fabs(qwt.value()) > SAF_bound ||
fabs(rhx.value()) > SAF_bound ||
fabs(rhy.value()) > SAF_bound ||
fabs(rhz.value()) > SAF_bound ||
fabs(rhw.value()) > SAF_bound ||
fabs(rwt.value()) > SAF_bound ||
fabs(shx.value()) > SAF_bound ||
fabs(shy.value()) > SAF_bound ||
fabs(shz.value()) > SAF_bound ||
fabs(shw.value()) > SAF_bound ||
fabs(swt.value()) > SAF_bound ||
fabs(thx.value()) > SAF_bound ||
fabs(thy.value()) > SAF_bound ||
fabs(thz.value()) > SAF_bound ||
fabs(thw.value()) > SAF_bound ||
fabs(twt.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(phx.value()));
SAF_bound = std::max(SAF_bound, fabs(phy.value()));
SAF_bound = std::max(SAF_bound, fabs(phz.value()));
SAF_bound = std::max(SAF_bound, fabs(phw.value()));
SAF_bound = std::max(SAF_bound, fabs(pwt.value()));
SAF_bound = std::max(SAF_bound, fabs(qhx.value()));
SAF_bound = std::max(SAF_bound, fabs(qhy.value()));
SAF_bound = std::max(SAF_bound, fabs(qhz.value()));
SAF_bound = std::max(SAF_bound, fabs(qhw.value()));
SAF_bound = std::max(SAF_bound, fabs(qwt.value()));
SAF_bound = std::max(SAF_bound, fabs(rhx.value()));
SAF_bound = std::max(SAF_bound, fabs(rhy.value()));
SAF_bound = std::max(SAF_bound, fabs(rhz.value()));
SAF_bound = std::max(SAF_bound, fabs(rhw.value()));
SAF_bound = std::max(SAF_bound, fabs(rwt.value()));
SAF_bound = std::max(SAF_bound, fabs(shx.value()));
SAF_bound = std::max(SAF_bound, fabs(shy.value()));
SAF_bound = std::max(SAF_bound, fabs(shz.value()));
SAF_bound = std::max(SAF_bound, fabs(shw.value()));
SAF_bound = std::max(SAF_bound, fabs(swt.value()));
SAF_bound = std::max(SAF_bound, fabs(thx.value()));
SAF_bound = std::max(SAF_bound, fabs(thy.value()));
SAF_bound = std::max(SAF_bound, fabs(thz.value()));
SAF_bound = std::max(SAF_bound, fabs(thw.value()));
SAF_bound = std::max(SAF_bound, fabs(twt.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) power_testH3_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return power_testH3_SAF(
Restricted_double(phx.value()),
Restricted_double(phy.value()),
Restricted_double(phz.value()),
Restricted_double(phw.value()),
Restricted_double(pwt.value()),
Restricted_double(qhx.value()),
Restricted_double(qhy.value()),
Restricted_double(qhz.value()),
Restricted_double(qhw.value()),
Restricted_double(qwt.value()),
Restricted_double(rhx.value()),
Restricted_double(rhy.value()),
Restricted_double(rhz.value()),
Restricted_double(rhw.value()),
Restricted_double(rwt.value()),
Restricted_double(shx.value()),
Restricted_double(shy.value()),
Restricted_double(shz.value()),
Restricted_double(shw.value()),
Restricted_double(swt.value()),
Restricted_double(thx.value()),
Restricted_double(thy.value()),
Restricted_double(thz.value()),
Restricted_double(thw.value()),
Restricted_double(twt.value()),
SAF_epsilon_0);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return power_testH3(
phx.exact(),
phy.exact(),
phz.exact(),
phw.exact(),
pwt.exact(),
qhx.exact(),
qhy.exact(),
qhz.exact(),
qhw.exact(),
qwt.exact(),
rhx.exact(),
rhy.exact(),
rhz.exact(),
rhw.exact(),
rwt.exact(),
shx.exact(),
shy.exact(),
shz.exact(),
shw.exact(),
swt.exact(),
thx.exact(),
thy.exact(),
thz.exact(),
thw.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -123,155 +451,6 @@ power_testH3(
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
/* */
Oriented_side
power_testH3(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &pwt,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qwt,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &rhx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &rhy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &rhz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &rhw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &rwt,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{
Oriented_side result = power_testH3(
phx.interval(),
phy.interval(),
phz.interval(),
phw.interval(),
pwt.interval(),
qhx.interval(),
qhy.interval(),
qhz.interval(),
qhw.interval(),
qwt.interval(),
rhx.interval(),
rhy.interval(),
rhz.interval(),
rhw.interval(),
rwt.interval(),
thx.interval(),
thy.interval(),
thz.interval(),
thw.interval(),
twt.interval());
FPU_set_cw(backup);
return result;
}
catch (Interval_nt_advanced::unsafe_comparison)
{
FPU_set_cw(backup);
return power_testH3(
phx.exact(),
phy.exact(),
phz.exact(),
phw.exact(),
pwt.exact(),
qhx.exact(),
qhy.exact(),
qhz.exact(),
qhw.exact(),
qwt.exact(),
rhx.exact(),
rhy.exact(),
rhz.exact(),
rhw.exact(),
rwt.exact(),
thx.exact(),
thy.exact(),
thz.exact(),
thw.exact(),
twt.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
/* */
Oriented_side
power_testH3(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &phw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &pwt,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qhw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &qwt,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thx,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thy,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thz,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &thw,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &twt)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{
Oriented_side result = power_testH3(
phx.interval(),
phy.interval(),
phz.interval(),
phw.interval(),
pwt.interval(),
qhx.interval(),
qhy.interval(),
qhz.interval(),
qhw.interval(),
qwt.interval(),
thx.interval(),
thy.interval(),
thz.interval(),
thw.interval(),
twt.interval());
FPU_set_cw(backup);
return result;
}
catch (Interval_nt_advanced::unsafe_comparison)
{
FPU_set_cw(backup);
return power_testH3(
phx.exact(),
phy.exact(),
phz.exact(),
phw.exact(),
pwt.exact(),
qhx.exact(),
qhy.exact(),
qhz.exact(),
qhw.exact(),
qwt.exact(),
thx.exact(),
thy.exact(),
thz.exact(),
thw.exact(),
twt.exact());
}
}
CGAL_END_NAMESPACE
#endif // CGAL_ARITHMETIC_FILTER_REGULAR_TRIANGULATION_RTH3_H

View File

@ -26,6 +26,152 @@
CGAL_BEGIN_NAMESPACE
inline
Oriented_side
in_smallest_orthogonalcircleC2_SAF(
const Static_filter_error &px,
const Static_filter_error &py,
const Static_filter_error &pw,
const Static_filter_error &qx,
const Static_filter_error &qy,
const Static_filter_error &qw,
const Static_filter_error &tx,
const Static_filter_error &ty,
const Static_filter_error &tw,
double & epsilon_0)
{
typedef Static_filter_error FT;
FT dpx = px-qx;
FT dpy = py-qy;
FT dtx = tx-qx;
FT dty = ty-qy;
FT dpz = square(dpx)+square(dpy);
return Oriented_side (sign_SAF((square(dtx)+square(dty)-tw+qw)*dpz
-(dpz-pw+qw)*(dpx*dtx+dpy*dty),
epsilon_0));
}
inline
Oriented_side
in_smallest_orthogonalcircleC2_SAF(
const Restricted_double &px,
const Restricted_double &py,
const Restricted_double &pw,
const Restricted_double &qx,
const Restricted_double &qy,
const Restricted_double &qw,
const Restricted_double &tx,
const Restricted_double &ty,
const Restricted_double &tw,
const double & epsilon_0)
{
typedef Restricted_double FT;
FT dpx = px-qx;
FT dpy = py-qy;
FT dtx = tx-qx;
FT dty = ty-qy;
FT dpz = square(dpx)+square(dpy);
return Oriented_side (sign_SAF((square(dtx)+square(dty)-tw+qw)*dpz
-(dpz-pw+qw)*(dpx*dtx+dpy*dty),
epsilon_0));
}
inline
Oriented_side
in_smallest_orthogonalcircleC2(
const Static_adaptatif_filter &px,
const Static_adaptatif_filter &py,
const Static_adaptatif_filter &pw,
const Static_adaptatif_filter &qx,
const Static_adaptatif_filter &qy,
const Static_adaptatif_filter &qw,
const Static_adaptatif_filter &tx,
const Static_adaptatif_filter &ty,
const Static_adaptatif_filter &tw)
{
bool re_adjusted = false;
static double SAF_bound = -1.0;
static double SAF_epsilon_0;
// Check the bounds. All arguments must be <= SAF_bound.
if ( fabs(px.value()) > SAF_bound ||
fabs(py.value()) > SAF_bound ||
fabs(pw.value()) > SAF_bound ||
fabs(qx.value()) > SAF_bound ||
fabs(qy.value()) > SAF_bound ||
fabs(qw.value()) > SAF_bound ||
fabs(tx.value()) > SAF_bound ||
fabs(ty.value()) > SAF_bound ||
fabs(tw.value()) > SAF_bound)
{
re_adjust:
// Re-adjust SAF_bound.
SAF_bound = std::max(0.0, fabs(px.value()));
SAF_bound = std::max(SAF_bound, fabs(py.value()));
SAF_bound = std::max(SAF_bound, fabs(pw.value()));
SAF_bound = std::max(SAF_bound, fabs(qx.value()));
SAF_bound = std::max(SAF_bound, fabs(qy.value()));
SAF_bound = std::max(SAF_bound, fabs(qw.value()));
SAF_bound = std::max(SAF_bound, fabs(tx.value()));
SAF_bound = std::max(SAF_bound, fabs(ty.value()));
SAF_bound = std::max(SAF_bound, fabs(tw.value()));
// recompute the epsilons: "just" call it over Static_filter_error.
// That's the tricky part that might not work for everything.
(void) in_smallest_orthogonalcircleC2_SAF(
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
Static_filter_error(SAF_bound),
SAF_epsilon_0);
// TODO: We should verify that all epsilons have really been updated.
}
try // Try the epsilon variant of the predicate.
{
return in_smallest_orthogonalcircleC2_SAF(
Restricted_double(px.value()),
Restricted_double(py.value()),
Restricted_double(pw.value()),
Restricted_double(qx.value()),
Restricted_double(qy.value()),
Restricted_double(qw.value()),
Restricted_double(tx.value()),
Restricted_double(ty.value()),
Restricted_double(tw.value()),
SAF_epsilon_0);
}
catch (Restricted_double::unsafe_comparison)
{
// It failed, we re-adjust once.
if (!re_adjusted) {
re_adjusted = true;
goto re_adjust;
}
// This scheme definitely fails => exact computation (filtered_exact<> ?).
return in_smallest_orthogonalcircleC2(
px.exact(),
py.exact(),
pw.exact(),
qx.exact(),
qy.exact(),
qw.exact(),
tx.exact(),
ty.exact(),
tw.exact());
}
}
#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION
template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE >
#endif
@ -42,6 +188,7 @@ in_smallest_orthogonalcircleC2(
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &ty,
const Filtered_exact <CGAL_IA_CT, CGAL_IA_ET, CGAL_IA_CACHE> &tw)
{
CGAL_assertion(Interval_nt_advanced::want_exceptions);
FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up);
try
{