added the first version for the files in the predicates subdir

This commit is contained in:
Menelaos Karavelas 2003-07-23 03:52:48 +00:00
parent a09d8538b2
commit 78db95fcb9
4 changed files with 3190 additions and 0 deletions

View File

@ -0,0 +1,725 @@
#ifndef CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H
#define CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H
#include <CGAL/predicates/Segment_Voronoi_diagram_predicates_C2.h>
#include <CGAL/predicates/check_filter.h>
CGAL_BEGIN_NAMESPACE
template<class FT, class Method_tag>
Comparison_result
svd_compare_distanceC2(const FT& qx, const FT& qy,
const FT& sx, const FT& sy,
const FT& tx, const FT& ty,
const FT& px, const FT& py, Method_tag)
{
// first check if (qx,qy) is inside, the boundary or at the exterior
// of the band of the segment s
must_be_filtered(qx);
FT dx = sx - tx;
FT dy = sy - ty;
FT d1x = qx - sx;
FT d1y = qy - sy;
FT d2x = qx - tx;
FT d2y = qy - ty;
if ( px == sx && py == sy ) {
FT o = dx * d1x + dy * d1y;
if ( o >= 0 ) {
return LARGER;
}
}
if ( px == tx && py == ty ) {
FT o = dx * d2x + dy * d2y;
if ( o <= 0 ) {
return LARGER;
}
}
FT d2_from_p = CGAL_NTS square (qx-px) + CGAL_NTS square(qy-py);
FT dot1 = dx * d1x + dy * d1y;
if ( dot1 >= 0 ) {
// q is outside (or the boundary of) the band on the side of s.
FT d2_from_s = CGAL_NTS square(d1x) + CGAL_NTS square(d1y);
return CGAL_NTS compare(d2_from_s, d2_from_p);
}
FT dot2 = dx * d2x + dy * d2y;
if ( dot2 <= 0 ) {
// q is outside (or the boundary of) the band on the side of t.
FT d2_from_t = CGAL_NTS square(d2x) + CGAL_NTS square(d2y);
return CGAL_NTS compare(d2_from_t, d2_from_p);
}
// q is strictly in the interior of the band, so I have to compare
// its distance from the supporting line.
FT c = sx * ty - sy * tx;
FT n = CGAL_NTS square(dx) + CGAL_NTS square(dy);
FT d2_from_l = CGAL_NTS square(qx * dy - qy * dx + c);
return CGAL_NTS compare(d2_from_l, d2_from_p * n);
}
template<class FT, class Method_tag>
Comparison_result
svd_compare_distanceC2(const FT& qx, const FT& qy,
const FT& s1x, const FT& s1y,
const FT& t1x, const FT& t1y,
const FT& s2x, const FT& s2y,
const FT& t2x, const FT& t2y, Method_tag)
{
// first check if (qx,qy) is inside, the boundary or at the exterior
// of the band of the segments s1, s2
must_be_filtered(qx);
FT d1x = s1x - t1x;
FT d1y = s1y - t1y;
FT d2x = s2x - t2x;
FT d2y = s2y - t2y;
FT dqs1x = qx - s1x;
FT dqs1y = qy - s1y;
FT dqt1x = qx - t1x;
FT dqt1y = qy - t1y;
FT dqs2x = qx - s2x;
FT dqs2y = qy - s2y;
FT dqt2x = qx - t2x;
FT dqt2y = qy - t2y;
FT dot_s1 = d1x * dqs1x + d1y * dqs1y;
int idx1; // 0 for s1, 1 for interior of 1, 2 for t1;
if ( qx == s1x && qy == s1y ) {
idx1 = 0;
} else if ( qx == t1x && qy == t1y ) {
idx1 = 2;
} else if ( dot_s1 >= 0 ) {
// q is outside (or the boundary of) the band of 1 on the side of s1.
idx1 = 0;
} else {
FT dot_t1 = d1x * dqt1x + d1y * dqt1y;
if ( dot_t1 <= 0 ) {
// q is outside (or the boundary of) the band of 1 on the side of t1.
idx1 = 2;
} else {
idx1 = 1;
}
}
FT dot_s2 = d2x * dqs2x + d2y * dqs2y;
int idx2; // 0 for s2, 1 for interior of 2, 2 for t2;
if ( qx == s2x && qy == s2y ) {
idx2 = 0;
} else if ( qx == t2x && qy == t2y ) {
idx2 = 2;
} else if ( dot_s2 >= 0 ) {
// q is outside (or the boundary of) the band of 2 on the side of s2.
idx2 = 0;
} else {
FT dot_t2 = d2x * dqt2x + d2y * dqt2y;
if ( dot_t2 <= 0 ) {
// q is outside (or the boundary of) the band of 2 on the side of t2.
idx2 = 2;
} else {
idx2 = 1;
}
}
if ( idx1 == 0 ) {
FT d2_from_s1 = CGAL_NTS square(dqs1x) + CGAL_NTS square(dqs1y);
// if ( qx == s1x && qy == s1y ) { d2_from_s1 = FT(0); }
if ( idx2 == 0 ) {
FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y);
// if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); }
if ( s1x == s2x && s1y == s2y ) { return EQUAL; }
// std::cout << "checkpoint 2.1" << std::endl;
return CGAL_NTS compare(d2_from_s1, d2_from_s2);
} else if ( idx2 == 2 ) {
// std::cout << "checkpoint 2.2" << std::endl;
FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y);
// if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); }
if ( s1x == t2x && s1y == t2y ) { return EQUAL; }
return CGAL_NTS compare(d2_from_s1, d2_from_t2);
} else { // idx2 == 1
FT c2 = s2x * t2y - s2y * t2x;
FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y);
FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2);
// std::cout << "checkpoint 2.3" << std::endl;
return CGAL_NTS compare(d2_from_s1 * n2, d2_from_l2);
}
} else if ( idx1 == 2 ) {
FT d2_from_t1 = CGAL_NTS square(dqt1x) + CGAL_NTS square(dqt1y);
// if ( qx == t1x && qy == t1y ) { d2_from_t1 = FT(0); }
if ( idx2 == 0 ) {
FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y);
// if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); }
/*
std::cout << "checkpoint 3.1" << std::endl;
std::pair<double,double> ivl;
ivl = to_interval(d2_from_t1);
std::cout << "interval of d2_from_t1: "
<< ivl.first << " " << ivl.second << std::endl;
ivl = to_interval(d2_from_s2);
std::cout << "interval of d2_from_s2: "
<< ivl.first << " " << ivl.second << std::endl;
*/
if ( t1x == s2x && t1y == s2y ) { return EQUAL; }
return CGAL_NTS compare(d2_from_t1, d2_from_s2);
} else if ( idx2 == 2 ) {
FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y);
// if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); }
// std::cout << "checkpoint 3.2" << std::endl;
if ( t1x == t2x && t1y == t2y ) { return EQUAL; }
return CGAL_NTS compare(d2_from_t1, d2_from_t2);
} else { // idx2 == 1
FT c2 = s2x * t2y - s2y * t2x;
FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y);
FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2);
// std::cout << "checkpoint 3.3" << std::endl;
return CGAL_NTS compare(d2_from_t1 * n2, d2_from_l2);
}
} else { // idx1 == 1
FT c1 = s1x * t1y - s1y * t1x;
FT n1 = CGAL_NTS square(d1x) + CGAL_NTS square(d1y);
FT d2_from_l1 = CGAL_NTS square(qx * d1y - qy * d1x + c1);
if ( idx2 == 0 ) {
FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y);
// if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); }
// std::cout << "checkpoint 4.1" << std::endl;
return CGAL_NTS compare(d2_from_l1, d2_from_s2 * n1);
} else if ( idx2 == 2 ) {
FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y);
// if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); }
// std::cout << "checkpoint 4.2" << std::endl;
return CGAL_NTS compare(d2_from_l1, d2_from_t2 * n1);
} else { // idx2 == 1
FT c2 = s2x * t2y - s2y * t2x;
FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y);
FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2);
// std::cout << "checkpoint 4.3" << std::endl;
return CGAL_NTS compare(d2_from_l1 * n2, d2_from_l2 * n1);
}
}
}
//--------------------------------------------------------------------------
#if 0
template<class FT, class Method_tag>
inline
Edge_conflict_test_result
svd_edge_predicate_ftC2(const std::vector<FT>& v,
char site_types[], int num_sites,
bool is_degenerate, Method_tag)
{
FT a;
must_be_filtered(a);
typedef Simple_cartesian<FT> Rep;
typedef Point_2<Rep> Point;
typedef Segment_2<Rep> Segment;
typedef Segment_Voronoi_diagram_site_2<Rep> Site;
Site t[num_sites];
for (int i = 0, k = 0; i < num_sites; i++) {
if ( site_types[i] == 'p' ) {
Point p(v[k], v[k+1]);
t[i].set_point( p );
} else if ( site_types[i] == 's' ) {
Point p1(v[k], v[k+1]), p2(v[k+2], v[k+3]);
Segment s(p1, p2);
t[i].set_segment( s );
} else {
CGAL_assertion( site_types[i] == 'p' ||
site_types[i] == 's' );
}
k += ( (site_types[i] == 'p') ? 2 : 4 );
}
CGAL_precondition( num_sites >= 3 && num_sites <= 5 );
if ( num_sites == 3 ) {
return svd_edge_test_degenerated_2(t[0], t[1], t[2]);
}
if ( num_sites == 4 ) {
if ( is_degenerate ) {
return svd_edge_test_degenerated_2(t[0], t[1], t[2], t[3]);
} else {
return svd_infinite_edge_test_2(t[0], t[1], t[2], t[3]);
}
}
return svd_edge_test_2(t[0], t[1], t[2], t[3], t[4]);
}
#endif
//--------------------------------------------------------------------------
template<class K, class Method_tag>
inline
Sign
svd_vertex_conflict_ftC2(const typename K::Site_2 t[],
int num_sites, Method_tag mtag)
{
typedef typename K::FT FT;
char site_types[num_sites];
std::vector<FT> v;
for (int i = 0; i < num_sites; i++) {
if ( t[i].is_point() ) {
v.push_back( t[i].point().x() );
v.push_back( t[i].point().y() );
site_types[i] = 'p';
} else {
v.push_back( t[i].source().x() );
v.push_back( t[i].source().y() );
v.push_back( t[i].target().x() );
v.push_back( t[i].target().y() );
site_types[i] = 's';
}
}
return svd_vertex_conflict_ftC2(v, site_types, num_sites, mtag);
}
template<class K, class Method_tag>
inline
Sign
svd_vertex_conflict_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& t,
Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, t};
return
svd_vertex_conflict_ftC2<K,Method_tag>(site_vec, 3, mtag);
}
template<class K, class Method_tag>
inline
Sign
svd_vertex_conflict_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& r,
const typename K::Site_2& t,
Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, r, t};
return
svd_vertex_conflict_ftC2<K,Method_tag>(site_vec, 4, mtag);
}
template<class FT, class Method_tag>
inline
Sign
svd_vertex_conflict_ftC2(const std::vector<FT>& v,
char site_types[], int num_sites,
Method_tag)
{
CGAL_precondition( num_sites == 3 || num_sites == 4 );
must_be_filtered(FT());
typedef Simple_cartesian<FT> Rep;
typedef CGAL::Segment_Voronoi_diagram_kernel_wrapper_2<Rep> Kernel;
typedef typename Kernel::Point_2 Point_2;
typedef typename Kernel::Segment_2 Segment_2;
typedef typename Kernel::Site_2 Site_2;
typedef Svd_incircle_2<Kernel,Method_tag> Incircle;
Site_2 t[num_sites];
for (int i = 0, k = 0; i < num_sites; i++) {
if ( site_types[i] == 'p' ) {
Point_2 p(v[k], v[k+1]);
t[i].set_point( p );
} else if ( site_types[i] == 's' ) {
Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]);
Segment_2 s(p1, p2);
t[i].set_segment( s );
} else {
CGAL_assertion( site_types[i] == 'p' ||
site_types[i] == 's' );
}
k += ( (site_types[i] == 'p') ? 2 : 4 );
}
if ( num_sites == 3 ) {
return Incircle()(t[0], t[1], t[2]);
}
return Incircle()(t[0], t[1], t[2], t[3]);
}
//--------------------------------------------------------------------------
template<class K, class Method_tag>
inline
bool
svd_finite_edge_conflict_ftC2(const typename K::Site_2 t[],
Sign sgn, int num_sites, Method_tag mtag)
{
typedef typename K::FT FT;
char site_types[num_sites];
std::vector<FT> v;
for (int i = 0; i < num_sites; i++) {
if ( t[i].is_point() ) {
v.push_back( t[i].point().x() );
v.push_back( t[i].point().y() );
site_types[i] = 'p';
} else {
v.push_back( t[i].source().x() );
v.push_back( t[i].source().y() );
v.push_back( t[i].target().x() );
v.push_back( t[i].target().y() );
site_types[i] = 's';
}
}
return svd_finite_edge_conflict_ftC2(v, sgn, site_types,
num_sites, mtag);
}
template<class K, class Method_tag>
inline
bool
svd_finite_edge_conflict_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& t,
Sign sgn, Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, t};
return
svd_finite_edge_conflict_ftC2<K,Method_tag>(site_vec, sgn, 3, mtag);
}
template<class K, class Method_tag>
inline
bool
svd_finite_edge_conflict_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& r,
const typename K::Site_2& t,
Sign sgn, Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, r, t};
return
svd_finite_edge_conflict_ftC2<K,Method_tag>(site_vec, sgn, 4, mtag);
}
template<class K, class Method_tag>
inline
bool
svd_finite_edge_conflict_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& r,
const typename K::Site_2& s,
const typename K::Site_2& t,
Sign sgn, Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, r, s, t};
return
svd_finite_edge_conflict_ftC2<K,Method_tag>(site_vec, sgn, 5, mtag);
}
template<class FT, class Method_tag>
inline
bool
svd_finite_edge_conflict_ftC2(const std::vector<FT>& v, Sign sgn,
char site_types[], int num_sites,
Method_tag)
{
CGAL_precondition( num_sites >= 3 || num_sites <= 5 );
must_be_filtered(FT());
typedef Simple_cartesian<FT> Rep;
typedef Segment_Voronoi_diagram_kernel_wrapper_2<Rep> Kernel;
typedef typename Kernel::Point_2 Point_2;
typedef typename Kernel::Segment_2 Segment_2;
typedef typename Kernel::Site_2 Site_2;
typedef Svd_finite_edge_interior_2<Kernel,Method_tag> Edge_interior_test;
Site_2 t[num_sites];
for (int i = 0, k = 0; i < num_sites; i++) {
if ( site_types[i] == 'p' ) {
Point_2 p(v[k], v[k+1]);
t[i].set_point( p );
} else if ( site_types[i] == 's' ) {
Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]);
Segment_2 s(p1, p2);
t[i].set_segment( s );
} else {
CGAL_assertion( site_types[i] == 'p' ||
site_types[i] == 's' );
}
k += ( (site_types[i] == 'p') ? 2 : 4 );
}
if ( num_sites == 3 ) {
return Edge_interior_test()(t[0], t[1], t[2], sgn);
} else if ( num_sites == 4 ) {
return Edge_interior_test()(t[0], t[1], t[2], t[3], sgn);
}
return Edge_interior_test()(t[0], t[1], t[2], t[3], t[4], sgn);
}
//--------------------------------------------------------------------------
template<class K, class Method_tag>
inline
bool
svd_infinite_edge_conflict_ftC2(const typename K::Site_2 t[],
Sign sgn, int num_sites, Method_tag mtag)
{
typedef typename K::FT FT;
char site_types[num_sites];
std::vector<FT> v;
for (int i = 0; i < num_sites; i++) {
if ( t[i].is_point() ) {
v.push_back( t[i].point().x() );
v.push_back( t[i].point().y() );
site_types[i] = 'p';
} else {
v.push_back( t[i].source().x() );
v.push_back( t[i].source().y() );
v.push_back( t[i].target().x() );
v.push_back( t[i].target().y() );
site_types[i] = 's';
}
}
return svd_infinite_edge_conflict_ftC2(v, sgn, site_types,
num_sites, mtag);
}
template<class K, class Method_tag>
inline
bool
svd_infinite_edge_conflict_ftC2(const typename K::Site_2& q,
const typename K::Site_2& r,
const typename K::Site_2& s,
const typename K::Site_2& t,
Sign sgn, Method_tag mtag)
{
typename K::Site_2 site_vec[] = {q, r, s, t};
return
svd_infinite_edge_conflict_ftC2<K,Method_tag>(site_vec, sgn, 4, mtag);
}
template<class FT, class Method_tag>
inline
bool
svd_infinite_edge_conflict_ftC2(const std::vector<FT>& v, Sign sgn,
char site_types[], int num_sites,
Method_tag)
{
CGAL_precondition( num_sites == 4 );
must_be_filtered(FT());
typedef Simple_cartesian<FT> Rep;
typedef Segment_Voronoi_diagram_kernel_wrapper_2<Rep> Kernel;
typedef typename Kernel::Point_2 Point_2;
typedef typename Kernel::Segment_2 Segment_2;
typedef typename Kernel::Site_2 Site_2;
typedef
Svd_infinite_edge_interior_2<Kernel,Method_tag> Edge_interior_test;
Site_2 t[num_sites];
for (int i = 0, k = 0; i < num_sites; i++) {
if ( site_types[i] == 'p' ) {
Point_2 p(v[k], v[k+1]);
t[i].set_point( p );
} else if ( site_types[i] == 's' ) {
Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]);
Segment_2 s(p1, p2);
t[i].set_segment( s );
} else {
CGAL_assertion( site_types[i] == 'p' ||
site_types[i] == 's' );
}
k += ( (site_types[i] == 'p') ? 2 : 4 );
}
return Edge_interior_test()(t[0], t[1], t[2], t[3], sgn);
}
//--------------------------------------------------------------------------
template<class K, class Method_tag>
inline
bool
svd_is_degenerate_edge_ftC2(const typename K::Site_2 t[],
int num_sites, Method_tag mtag)
{
typedef typename K::FT FT;
char site_types[num_sites];
std::vector<FT> v;
for (int i = 0; i < num_sites; i++) {
if ( t[i].is_point() ) {
v.push_back( t[i].point().x() );
v.push_back( t[i].point().y() );
site_types[i] = 'p';
} else {
v.push_back( t[i].source().x() );
v.push_back( t[i].source().y() );
v.push_back( t[i].target().x() );
v.push_back( t[i].target().y() );
site_types[i] = 's';
}
}
return svd_is_degenerate_edge_ftC2(v, site_types,
num_sites, mtag);
}
template<class K, class Method_tag>
inline
bool
svd_is_degenerate_edge_ftC2(const typename K::Site_2& p,
const typename K::Site_2& q,
const typename K::Site_2& r,
const typename K::Site_2& t,
Method_tag mtag)
{
typename K::Site_2 site_vec[] = {p, q, r, t};
return svd_is_degenerate_edge_ftC2<K,Method_tag>(site_vec, 4, mtag);
}
template<class FT, class Method_tag>
inline
bool
svd_is_degenerate_edge_ftC2(const std::vector<FT>& v,
char site_types[], int num_sites,
Method_tag)
{
CGAL_precondition( num_sites == 4 );
must_be_filtered(FT());
typedef Simple_cartesian<FT> Rep;
typedef Segment_Voronoi_diagram_kernel_wrapper_2<Rep> Kernel;
typedef typename Kernel::Point_2 Point_2;
typedef typename Kernel::Segment_2 Segment_2;
typedef typename Kernel::Site_2 Site_2;
typedef
Svd_is_degenerate_edge_test_2<Kernel,Method_tag> Is_degenerate_edge;
Site_2 t[num_sites];
for (int i = 0, k = 0; i < num_sites; i++) {
if ( site_types[i] == 'p' ) {
Point_2 p(v[k], v[k+1]);
t[i].set_point( p );
} else if ( site_types[i] == 's' ) {
Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]);
Segment_2 s(p1, p2);
t[i].set_segment( s );
} else {
CGAL_assertion( site_types[i] == 'p' ||
site_types[i] == 's' );
}
k += ( (site_types[i] == 'p') ? 2 : 4 );
}
return Is_degenerate_edge()(t[0], t[1], t[2], t[3]);
}
//--------------------------------------------------------------------------
CGAL_END_NAMESPACE
#ifdef CGAL_ARITHMETIC_FILTER_H
#ifndef CGAL_ARITHMETIC_FILTER_SVD_PREDICATES_FTC2_H
#include <CGAL/Arithmetic_filter/predicates/svd_predicates_ftC2.h>
#endif // CGAL_ARITHMETIC_FILTER_SVD_PREDICATES_FTC2_H
#endif
#endif // CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H

View File

@ -0,0 +1,327 @@
#ifndef CGAL_SQUARE_ROOT_1_H
#define CGAL_SQUARE_ROOT_1_H
#include <CGAL/basic.h>
#include <CGAL/enum.h>
#include <iostream>
#define CHECK_CGAL_PRECONDITIONS 0
CGAL_BEGIN_NAMESPACE
template<class NT>
class Square_root_1
{
private:
NT x, y, r;
private:
typedef Square_root_1<NT> Self;
public:
typedef NT FT;
typedef NT RT;
Square_root_1() : x(0), y(0), r(0) {}
Square_root_1(int i) : x(i), y(0), r(0) {}
Square_root_1(const NT& a) : x(a), y(0), r(0) {}
Square_root_1(const NT& a, const NT& b, const NT& c)
: x(a), y(b), r(c)
{
// CGAL_assertion( !(CGAL_NTS is_negative(r)) );
}
Square_root_1(const Square_root_1<NT>& other)
: x(other.x), y(other.y), r(other.r) {}
NT a() const { return x; }
NT b() const { return y; }
NT c() const { return r; }
Self square() const
{
NT xy = x * y;
return Self(CGAL_NTS square(x) + CGAL_NTS square(y) * r,
xy + xy,
r);
}
Sign sign() const
{
Sign sx = CGAL_NTS sign(x);
if ( CGAL_NTS sign(r) == ZERO ) { return sx; }
Sign sy = CGAL_NTS sign(y);
if ( sx == sy ) { return sx; }
if ( sx == ZERO ) { return sy; }
return Sign( sx * CGAL_NTS compare( CGAL_NTS square(x),
r * CGAL_NTS square(y) )
);
}
Self operator+() const
{
return (*this);
}
Self operator-() const
{
return Self(-x, -y, r);
}
double to_double() const
{
// THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS
double xd = CGAL_NTS to_double(x);
double yd = CGAL_NTS to_double(y);
double rd = CGAL_NTS to_double(r);
return (xd + yd * CGAL_NTS sqrt(rd));
}
std::pair<double,double> to_interval() const
{
// THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS
std::pair<double,double> x_ivl = CGAL::to_interval(x);
std::pair<double,double> y_ivl = CGAL::to_interval(y);
std::pair<double,double> r_ivl = CGAL::to_interval(r);
std::pair<double,double> sqrt_r_ivl(CGAL_NTS sqrt(r_ivl.first),
CGAL_NTS sqrt(r_ivl.second));
std::pair<double,double>
ivl(x_ivl.first + y_ivl.first * sqrt_r_ivl.first,
x_ivl.second + y_ivl.second * sqrt_r_ivl.second);
return ivl;
}
};
// operator *
template<class NT>
inline
Square_root_1<NT>
operator*(const Square_root_1<NT>& x, const NT& n)
{
return Square_root_1<NT>(x.a() * n, x.b() * n, x.c());
}
template<class NT>
inline
Square_root_1<NT>
operator*(const NT& n, const Square_root_1<NT>& x)
{
return (x * n);
}
template<class NT>
inline
Square_root_1<NT>
operator*(const Square_root_1<NT>& x, const Square_root_1<NT>& y)
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL );
#endif
NT a = x.a() * y.a() + x.b() * y.b() * x.c();
NT b = x.a() * y.b() + x.b() * y.a();
return Square_root_1<NT>(a, b, x.c());
}
// operator +
template<class NT>
inline
Square_root_1<NT>
operator+(const Square_root_1<NT>& x, const NT& n)
{
return Square_root_1<NT>(x.a() + n, x.b(), x.c());
}
template<class NT>
inline
Square_root_1<NT>
operator+(const NT& n, const Square_root_1<NT>& x)
{
return (x + n);
}
template<class NT>
inline
Square_root_1<NT>
operator+(const Square_root_1<NT>& x, const Square_root_1<NT>& y)
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL );
#endif
return Square_root_1<NT>(x.a() + y.a(), x.b() + y.b(), x.c());
}
// operator -
template<class NT>
inline
Square_root_1<NT>
operator-(const Square_root_1<NT>& x, const NT& n)
{
return x + (-n);
}
template<class NT>
inline
Square_root_1<NT>
operator-(const NT& n, const Square_root_1<NT>& x)
{
return -(x - n);
}
template<class NT>
inline
Square_root_1<NT>
operator-(const Square_root_1<NT>& x, const Square_root_1<NT>& y)
{
return (x + (-y));
}
template<class NT>
inline
std::pair<double,double>
to_interval(const Square_root_1<NT>& x)
{
return x.to_interval();
}
namespace NTS {
template<class NT>
inline
bool
is_positive(const Square_root_1<NT>& x)
{
return sign(x) == POSITIVE;
}
template<class NT>
inline
bool
is_negative(const Square_root_1<NT>& x)
{
return sign(x) == NEGATIVE;
}
template<class NT>
inline
bool
is_zero(const Square_root_1<NT>& x)
{
return sign(x) == ZERO;
}
template<class NT>
inline
Sign
sign(const Square_root_1<NT>& x)
{
return x.sign();
}
template<class NT>
inline
Square_root_1<NT>
square(const Square_root_1<NT>& x)
{
return x.square();
}
template<class NT>
inline
Comparison_result
compare(const Square_root_1<NT>& x, const Square_root_1<NT>& y)
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL );
#endif
Sign s = CGAL_NTS sign(x - y);
if ( s == ZERO ) { return EQUAL; }
return (s == POSITIVE) ? LARGER : SMALLER;
}
template<class NT>
inline
double
to_double(const Square_root_1<NT>& x)
{
return x.to_double();
}
} // namespace NTS
// operator <<
template<class Stream, class NT>
inline
Stream&
operator<<(Stream& os, const Square_root_1<NT>& x)
{
os << "(" << x.a() << ")+(" << x.b() << ") sqrt{" << x.c() << "}";
return os;
#if 0
if ( CGAL_NTS is_zero(x) ) {
os << "0";
return os;
}
// Square_root_1<NT> One(NT(1), NT(0), x.r());
if ( CGAL_NTS sign(x.a()) != ZERO ) {
os << x.a();
if ( CGAL_NTS is_positive(x.b()) ) {
os << "+";
}
}
if ( CGAL_NTS sign(x.b()) != ZERO &&
CGAL_NTS sign(x.c()) != ZERO ) {
// if ( CGAL_NTS sign(x.b() - One) != ZERO ) {
os << x.b() << " ";
// }
os << "sqrt{" << x.c() << "}";
}
return os;
#endif
}
CGAL_END_NAMESPACE
#endif // CGAL_SQUARE_ROOT_1_H

View File

@ -0,0 +1,379 @@
#ifndef CGAL_SQUARE_ROOT_2_H
#define CGAL_SQUARE_ROOT_2_H
#include <CGAL/predicates/Square_root_1.h>
CGAL_BEGIN_NAMESPACE
template<class NT>
class Square_root_2
{
private:
typedef Square_root_2<NT> Self;
typedef Square_root_1<NT> Sqrt_1;
NT a0_, a1_, a2_, a3_;
NT A_, B_;
public:
typedef NT FT;
typedef NT RT;
public:
Square_root_2()
: a0_(0), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {}
Square_root_2(int i)
: a0_(i), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {}
Square_root_2(const NT& a)
: a0_(a), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {}
Square_root_2(const NT& a0, const NT& a1, const NT& a2,
const NT& a3, const NT& A, const NT& B)
: a0_(a0), a1_(a1), a2_(a2), a3_(a3), A_(A), B_(B)
{
// CGAL_assertion( !(CGAL_NTS is_negative(A_)) );
// CGAL_assertion( !(CGAL_NTS is_negative(B_)) );
}
Square_root_2(const Square_root_2<NT>& other)
: a0_(other.a0_), a1_(other.a1_), a2_(other.a2_),
a3_(other.a3_), A_(other.A_), B_(other.B_) {}
NT a() const { return a0_; }
NT b() const { return a1_; }
NT c() const { return a2_; }
NT d() const { return a3_; }
NT e() const { return A_; }
NT f() const { return B_; }
Self operator*(const Self& b) const
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(A_, b.A_) == EQUAL );
CGAL_precondition( CGAL_NTS compare(B_, b.B_) == EQUAL );
#endif
NT a0 = a0_ * b.a0_ + a1_ * b.a1_ * A_ + a2_ * b.a2_ * B_
+ a3_ * b.a3_ * A_ * B_;
NT a1 = a0_ * b.a1_ + a1_ * b.a0_ + (a2_ * b.a3_ + a3_ * b.a2_) * B_;
NT a2 = a0_ * b.a2_ + a2_ * b.a0_ + (a1_ * b.a3_ + a3_ * b.a1_) * A_;
NT a3 = a0_ * b.a3_ + a3_ * b.a0_ + a1_ * b.a2_ + a2_ * b.a1_;
return Self(a0, a1, a2, a3, A_, B_);
}
Self operator-() const
{
return Self(-a0_, -a1_, -a2_, -a3_, A_, B_);
}
Self operator+() const
{
return (*this);
}
Self square() const
{
NT a0 = CGAL_NTS square(a0_) + CGAL_NTS square(a1_) * A_
+ CGAL_NTS square(a2_) * B_ + CGAL_NTS square(a3_) * A_ * B_;
NT a1_half = a0_ * a1_ + a2_ * a3_ * B_;
NT a2_half = a0_ * a2_ + a1_ * a3_ * A_;
NT a3_half = a0_ * a3_ + a1_ * a2_;
NT a1 = a1_half + a1_half;
NT a2 = a2_half + a2_half;
NT a3 = a3_half + a3_half;
return Self(a0, a1, a2, a3, A_, B_);
}
Sign sign() const
{
// std::cout << "cp sgn 1" << std::flush;
Sqrt_1 x(a0_, a1_, A_);
Sqrt_1 y(a2_, a3_, A_);
#if 0
std::cout << "a0: " << a0_ << std::endl;
std::cout << "a1: " << a1_ << std::endl;
std::cout << "A : " << A_ << std::endl;
std::cout << "x: " << x << std::endl;
std::cout << " 2" << std::flush;
#endif
Sign s_x = CGAL_NTS sign(x);
// std::cout << " 3" << std::flush;
Sign s_y = CGAL_NTS sign(y);
// std::cout << " 4" << std::flush;
Sign s_B = CGAL_NTS sign(B_);
// std::cout << " 5" << std::flush;
// std::cout << std::endl;
if ( s_B == ZERO ) {
return s_x;
} else if ( s_x == s_y ) {
return s_x;
} else if ( s_x == ZERO ) {
return s_y;
} else if ( s_y == ZERO ) {
return s_x;
} else {
Sqrt_1 Q = CGAL_NTS square(x) - CGAL_NTS square(y) * B_;
return Sign(s_x * CGAL_NTS sign(Q));
}
}
double to_double() const
{
// THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS
double a0d = CGAL_NTS to_double(a0_);
double a1d = CGAL_NTS to_double(a1_);
double a2d = CGAL_NTS to_double(a2_);
double a3d = CGAL_NTS to_double(a3_);
double Ad = CGAL_NTS to_double(A_);
double Bd = CGAL_NTS to_double(B_);
return (a0d + a1d * CGAL_NTS sqrt(Ad) + a2d * CGAL_NTS sqrt(Bd)
+ a3d * CGAL_NTS sqrt(Ad * Bd));
}
};
// operator *
template<class NT>
inline
Square_root_2<NT>
operator*(const Square_root_2<NT>& x, const NT& n)
{
return Square_root_2<NT>(x.a() * n, x.b() * n, x.c() * n,
x.d() * n, x.e(), x.f());
}
template<class NT>
inline
Square_root_2<NT>
operator*(const NT& n, const Square_root_2<NT>& x)
{
return (x * n);
}
// operator +
template<class NT>
inline
Square_root_2<NT>
operator+(const Square_root_2<NT>& x, const NT& n)
{
return Square_root_2<NT>(x.a() + n, x.b(), x.c(), x.d(),
x.e(), x.f());
}
template<class NT>
inline
Square_root_2<NT>
operator+(const NT& n, const Square_root_2<NT>& x)
{
return (x + n);
}
template<class NT>
inline
Square_root_2<NT>
operator+(const Square_root_2<NT>& x, const Square_root_2<NT>& y)
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(x.e(), y.e()) == EQUAL );
CGAL_precondition( CGAL_NTS compare(x.f(), y.f()) == EQUAL );
#endif
return Square_root_2<NT>(x.a() + y.a(), x.b() + y.b(),
x.c() + y.c(), x.d() + y.d(),
x.e(), x.f());
}
// operator -
template<class NT>
inline
Square_root_2<NT>
operator-(const Square_root_2<NT>& x, const NT& n)
{
return x + (-n);
}
template<class NT>
inline
Square_root_2<NT>
operator-(const NT& n, const Square_root_2<NT>& x)
{
return -(x - n);
}
template<class NT>
inline
Square_root_2<NT>
operator-(const Square_root_2<NT>& x, const Square_root_2<NT>& y)
{
return (x + (-y));
}
namespace NTS {
template<class NT>
inline
bool
is_positive(const Square_root_2<NT>& x)
{
return sign(x) == POSITIVE;
}
template<class NT>
inline
bool
is_negative(const Square_root_2<NT>& x)
{
return sign(x) == NEGATIVE;
}
template<class NT>
inline
bool
is_zero(const Square_root_2<NT>& x)
{
return sign(x) == ZERO;
}
template<class NT>
inline
Sign
sign(const Square_root_2<NT>& x)
{
return x.sign();
}
template<class NT>
inline
Square_root_2<NT>
square(const Square_root_2<NT>& x)
{
return x.square();
}
template<class NT>
inline
Comparison_result
compare(const Square_root_2<NT>& x,
const Square_root_2<NT>& y)
{
#if CHECK_CGAL_PRECONDITIONS
CGAL_precondition( CGAL_NTS compare(x.e(), y.e()) == EQUAL );
CGAL_precondition( CGAL_NTS compare(x.f(), y.f()) == EQUAL );
#endif
Sign s = CGAL_NTS sign(x - y);
if ( s == ZERO ) { return EQUAL; }
return (s == POSITIVE) ? LARGER : SMALLER;
}
} // namespace NTS
// operator <<
template<class Stream, class NT>
inline
Stream&
operator<<(Stream& os, const Square_root_2<NT>& x)
{
os << "(" << x.a() << ")+(" << x.b() << ") sqrt{" << x.e() << "}";
os << "+(" << x.c() << ") sqrt{" << x.f() << "}";
os << "+(" << x.d() << ") sqrt{(" << x.e() << ") (" << x.f()
<< ")}";
return os;
#if 0
if ( CGAL_NTS is_zero(x) ) {
os << "0";
return os;
}
Square_root_2<NT> One(NT(1), NT(0), NT(0), NT(0), x.e(), x.f());
if ( CGAL_NTS sign(x.a()) != ZERO ) {
os << x.a();
if ( CGAL_NTS is_positive(x.b()) ) {
os << "+";
}
}
if ( CGAL_NTS sign(x.b()) != ZERO &&
CGAL_NTS sign(x.e()) != ZERO ) {
if ( CGAL_NTS sign(x.b() - One) != ZERO ) {
os << x.b() << " ";
}
os << "sqrt{" << x.e() << "}";
if ( CGAL_NTS is_positive(x.c()) ) {
os << "+";
}
}
if ( CGAL_NTS sign(x.c()) != ZERO &&
CGAL_NTS sign(x.f()) != ZERO ) {
if ( CGAL_NTS sign(x.c() - One) != ZERO ) {
os << x.c() << " ";
}
os << "sqrt{" << x.f() << "}";
if ( CGAL_NTS is_positive(x.d()) ) {
os << "+";
}
}
if ( CGAL_NTS sign(x.d()) != ZERO &&
CGAL_NTS sign(x.e()) != ZERO &&
CGAL_NTS sign(x.f()) != ZERO ) {
if ( CGAL_NTS sign(x.d() - One) != ZERO ) {
os << x.d() << " ";
}
os << "sqrt{" << x.e() << " " << x.f() << "}";
}
return os;
#endif
}
CGAL_END_NAMESPACE
#endif // CGAL_SQUARE_ROOT_2_H