Cleaned up

This commit is contained in:
Efi Fogel 2022-05-18 13:37:17 +03:00
parent 418e2b0c6b
commit d1b05fd144
5 changed files with 683 additions and 850 deletions

View File

@ -17,4 +17,5 @@ endforeach()
if(CGAL_Qt5_FOUND) if(CGAL_Qt5_FOUND)
target_link_libraries(draw_arr PUBLIC CGAL::CGAL_Basic_viewer) target_link_libraries(draw_arr PUBLIC CGAL::CGAL_Basic_viewer)
target_link_libraries(conics PUBLIC CGAL::CGAL_Basic_viewer)
endif() endif()

View File

@ -2,6 +2,7 @@
// Constructing an arrangement of various conic arcs. // Constructing an arrangement of various conic arcs.
#include <CGAL/config.h> #include <CGAL/config.h>
#include <CGAL/draw_arrangement_2.h>
#ifdef CGAL_USE_CORE #ifdef CGAL_USE_CORE
@ -11,6 +12,7 @@
int main() { int main() {
Arrangement arr; Arrangement arr;
#if 0
// Insert a hyperbolic arc (C1), supported by the hyperbola y = 1/x // Insert a hyperbolic arc (C1), supported by the hyperbola y = 1/x
// (or: xy - 1 = 0) with the endpoints (1/4, 4) and (2, 1/2). // (or: xy - 1 = 0) with the endpoints (1/4, 4) and (2, 1/2).
// The arc is counterclockwise oriented. // The arc is counterclockwise oriented.
@ -20,21 +22,24 @@ int main() {
// Insert a full ellipse (C2), which is (x/4)^2 + (y/2)^2 = 0 rotated by // Insert a full ellipse (C2), which is (x/4)^2 + (y/2)^2 = 0 rotated by
// phi = 36.87 degrees (such that sin(phi) = 0.6, cos(phi) = 0.8), // phi = 36.87 degrees (such that sin(phi) = 0.6, cos(phi) = 0.8),
// yielding: 58x^2 + 72y^2 - 48xy - 360 = 0. // yielding: 58x^2 + 72y^2 - 48xy - 360 = 0.
insert(arr, Conic_arc (58, 72, -48, 0, 0, -360)); insert(arr, Conic_arc(58, 72, -48, 0, 0, -360));
#endif
// Insert the segment (C3) (1, 1) -- (0, -3). // Insert the segment (C3) (1, 1) -- (0, -3).
insert(arr, Conic_arc(Rat_segment(Rat_point(1, 1), Rat_point(0, -3)))); insert(arr, Conic_arc(Rat_segment(Rat_point(1, 1), Rat_point(0, -3))));
#if 0
// Insert a circular arc (C4) supported by the circle x^2 + y^2 = 5^2, // Insert a circular arc (C4) supported by the circle x^2 + y^2 = 5^2,
// with (-3, 4) and (4, 3) as its endpoints. We want the arc to be // with (-3, 4) and (4, 3) as its endpoints. We want the arc to be
// clockwise-oriented, so it passes through (0, 5) as well. // clockwise-oriented, so it passes through (0, 5) as well.
Conic_arc c4(Rat_point(-3, 4), Rat_point(0, 5), Rat_point(4, 3)); Conic_arc c4(Rat_point(-3, 4), Rat_point(0, 5), Rat_point(4, 3));
insert(arr, c4); insert(arr, c4);
#endif
// Insert a full unit circle (C5) that is centered at (0, 4). // Insert a full unit circle (C5) that is centered at (0, 4).
insert(arr, Conic_arc(Rat_circle(Rat_point(0,4), 1))); insert(arr, Conic_arc(Rat_circle(Rat_point(0,4), 1)));
#if 0
// Insert a parabolic arc (C6) supported by the parabola y = -x^2 with // Insert a parabolic arc (C6) supported by the parabola y = -x^2 with
// endpoints (-sqrt(3),-3) (~(-1.73,-3)) and (sqrt(2),-2) (~(1.41,-2)). // endpoints (-sqrt(3),-3) (~(-1.73,-3)) and (sqrt(2),-2) (~(1.41,-2)).
// Since the x-coordinates of the endpoints cannot be acccurately represented, // Since the x-coordinates of the endpoints cannot be acccurately represented,
@ -53,8 +58,12 @@ int main() {
// is 1/2 (therefore its squared radius is 1/4) (C7). // is 1/2 (therefore its squared radius is 1/4) (C7).
Rat_circle circ7(Rat_point(4, Rational(5,2)), Rational(1,4)); Rat_circle circ7(Rat_point(4, Rational(5,2)), Rational(1,4));
insert(arr, Conic_arc(circ7, CGAL::CLOCKWISE, Point(4, 3), Point(4, 2))); insert(arr, Conic_arc(circ7, CGAL::CLOCKWISE, Point(4, 3), Point(4, 2)));
#endif
print_arrangement_size(arr); print_arrangement_size(arr);
CGAL::draw(arr);
return 0; return 0;
} }
@ -62,8 +71,7 @@ int main() {
#include <iostream> #include <iostream>
int main () int main() {
{
std::cout << "Sorry, this example needs GMP and CORE\n"; std::cout << "Sorry, this example needs GMP and CORE\n";
return 0; return 0;
} }

View File

@ -628,31 +628,31 @@ public:
// Compute the x- and y-coordinates of intersection points of the base // Compute the x- and y-coordinates of intersection points of the base
// conic and the k'th auxiliary conic. // conic and the k'th auxiliary conic.
n_xs = _compute_resultant_roots(nt_traits, n_xs = compute_resultant_roots(nt_traits,
base_coeffs[0], base_coeffs[1], base_coeffs[0], base_coeffs[1],
base_coeffs[2], base_coeffs[2],
base_coeffs[3], base_coeffs[4], base_coeffs[3], base_coeffs[4],
base_coeffs[5], base_coeffs[5],
deg_base, deg_base,
aux_coeffs[0], aux_coeffs[1], aux_coeffs[0], aux_coeffs[1],
aux_coeffs[2], aux_coeffs[2],
aux_coeffs[3], aux_coeffs[4], aux_coeffs[3], aux_coeffs[4],
aux_coeffs[5], aux_coeffs[5],
deg_aux, deg_aux,
xs); xs);
n_ys = _compute_resultant_roots(nt_traits, n_ys = compute_resultant_roots(nt_traits,
base_coeffs[1], base_coeffs[0], base_coeffs[1], base_coeffs[0],
base_coeffs[2], base_coeffs[2],
base_coeffs[4], base_coeffs[3], base_coeffs[4], base_coeffs[3],
base_coeffs[5], base_coeffs[5],
deg_base, deg_base,
aux_coeffs[1], aux_coeffs[0], aux_coeffs[1], aux_coeffs[0],
aux_coeffs[2], aux_coeffs[2],
aux_coeffs[4], aux_coeffs[3], aux_coeffs[4], aux_coeffs[3],
aux_coeffs[5], aux_coeffs[5],
deg_aux, deg_aux,
ys); ys);
// Find the intersection point which is nearest the given approximation // Find the intersection point which is nearest the given approximation
// and set it as the endpoint. // and set it as the endpoint.
@ -667,8 +667,7 @@ public:
nt_traits.convert(base_coeffs[4]) * ys[j] + nt_traits.convert(base_coeffs[4]) * ys[j] +
nt_traits.convert(base_coeffs[5]); nt_traits.convert(base_coeffs[5]);
if (CGAL::sign(val) != ZERO) if (CGAL::sign(val) != ZERO) continue;
continue;
val = nt_traits.convert(aux_coeffs[0]) * xs[i]*xs[i] + val = nt_traits.convert(aux_coeffs[0]) * xs[i]*xs[i] +
nt_traits.convert(aux_coeffs[1]) * ys[j]*ys[j] + nt_traits.convert(aux_coeffs[1]) * ys[j]*ys[j] +
@ -749,7 +748,7 @@ public:
// Duplicate the extra data, if necessary. // Duplicate the extra data, if necessary.
m_extra_data = (arc.m_extra_data != nullptr) ? m_extra_data = (arc.m_extra_data != nullptr) ?
new Extra_data(*(arc.m_extra_data)) ? nullptr; new Extra_data(*(arc.m_extra_data)) : nullptr;
return (*this); return (*this);
} }
@ -764,12 +763,12 @@ public:
/*! Get the coefficients of the underlying conic. /*! Get the coefficients of the underlying conic.
*/ */
const Integer& r() const {return (m_r);} const Integer& r() const { return (m_r); }
const Integer& s() const {return (m_s);} const Integer& s() const { return (m_s); }
const Integer& t() const {return (m_t);} const Integer& t() const { return (m_t); }
const Integer& u() const {return (m_u);} const Integer& u() const { return (m_u); }
const Integer& v() const {return (m_v);} const Integer& v() const { return (m_v); }
const Integer& w() const {return (m_w);} const Integer& w() const { return (m_w); }
/*! Check whether the arc is x-monotone. /*! Check whether the arc is x-monotone.
*/ */

View File

@ -8,7 +8,7 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
// //
// //
// Author(s) : Ron Wein <wein@post.tau.ac.il> // Author(s): Ron Wein <wein@post.tau.ac.il>
#ifndef CGAL_CONIC_INTERSECTIONS_2_H #ifndef CGAL_CONIC_INTERSECTIONS_2_H
#define CGAL_CONIC_INTERSECTIONS_2_H #define CGAL_CONIC_INTERSECTIONS_2_H
@ -24,8 +24,7 @@
namespace CGAL { namespace CGAL {
/*! /*! Compute the roots of the resultants of the two bivariate polynomials:
* Compute the roots of the resultants of the two bivariate polynomials:
* C1: r1*x^2 + s1*y^2 + t1*xy + u1*x + v1*y + w1 = 0 * C1: r1*x^2 + s1*y^2 + t1*xy + u1*x + v1*y + w1 = 0
* C2: r2*x^2 + s2*y^2 + t2*xy + u2*x + v2*y + w2 = 0 * C2: r2*x^2 + s2*y^2 + t2*xy + u2*x + v2*y + w2 = 0
* \param deg1 The degree of the first curve. * \param deg1 The degree of the first curve.
@ -35,84 +34,76 @@ namespace CGAL {
* \pre xs must be a vector of size 4. * \pre xs must be a vector of size 4.
* \return The number of distinct roots found. * \return The number of distinct roots found.
*/ */
template <class Nt_traits> template <typename Nt_traits>
int int compute_resultant_roots(Nt_traits& nt_traits,
_compute_resultant_roots (Nt_traits& nt_traits, const typename Nt_traits::Integer& r1,
const typename Nt_traits::Integer& r1, const typename Nt_traits::Integer& s1,
const typename Nt_traits::Integer& s1, const typename Nt_traits::Integer& t1,
const typename Nt_traits::Integer& t1, const typename Nt_traits::Integer& u1,
const typename Nt_traits::Integer& u1, const typename Nt_traits::Integer& v1,
const typename Nt_traits::Integer& v1, const typename Nt_traits::Integer& w1,
const typename Nt_traits::Integer& w1, const int& deg1,
const int& deg1, const typename Nt_traits::Integer& r2,
const typename Nt_traits::Integer& r2, const typename Nt_traits::Integer& s2,
const typename Nt_traits::Integer& s2, const typename Nt_traits::Integer& t2,
const typename Nt_traits::Integer& t2, const typename Nt_traits::Integer& u2,
const typename Nt_traits::Integer& u2, const typename Nt_traits::Integer& v2,
const typename Nt_traits::Integer& v2, const typename Nt_traits::Integer& w2,
const typename Nt_traits::Integer& w2, const int& deg2,
const int& deg2, typename Nt_traits::Algebraic* xs)
typename Nt_traits::Algebraic *xs)
{ {
if (deg1 == 2 && deg2 == 1) if ((deg1 == 2) && (deg2 == 1)) {
{
// If necessary, swap roles between the two curves, so that the first // If necessary, swap roles between the two curves, so that the first
// curve always has the minimal degree. // curve always has the minimal degree.
return (_compute_resultant_roots (nt_traits, return (compute_resultant_roots(nt_traits,
r2, s2, t2, u2, v2, w2, r2, s2, t2, u2, v2, w2,
deg2, deg2,
r1, s1, t1, u1, v1, w1, r1, s1, t1, u1, v1, w1,
deg1, deg1,
xs)); xs));
} }
// Act according to the degree of the first conic curve. // Act according to the degree of the first conic curve.
const typename Nt_traits::Integer _two = 2; const typename Nt_traits::Integer two = 2;
typename Nt_traits::Integer c[5]; typename Nt_traits::Integer c[5];
unsigned int degree = 4; unsigned int degree = 4;
typename Nt_traits::Algebraic *xs_end; typename Nt_traits::Algebraic* xs_end;
if (deg1 == 1) if (deg1 == 1) {
{
// The first curve has no quadratic coefficients, and represents a line. // The first curve has no quadratic coefficients, and represents a line.
if (CGAL::sign (v1) == ZERO) if (CGAL::sign (v1) == ZERO) {
{
// The first line is u1*x + w1 = 0, therefore: // The first line is u1*x + w1 = 0, therefore:
xs[0] = nt_traits.convert(-w1) / nt_traits.convert(u1); xs[0] = nt_traits.convert(-w1) / nt_traits.convert(u1);
return (1); return 1;
} }
// We can write the first curve as: y = -(u1*x + w1) / v1. // We can write the first curve as: y = -(u1*x + w1) / v1.
if (deg2 == 1) if (deg2 == 1) {
{
// The second curve is also a line. We therefore get the linear // The second curve is also a line. We therefore get the linear
// equation c[1]*x + c[0] = 0: // equation c[1]*x + c[0] = 0:
c[1] = v1*u2 - u1*v2; c[1] = v1*u2 - u1*v2;
c[0] = v1*w2 - w1*v2; c[0] = v1*w2 - w1*v2;
if (CGAL::sign (c[1]) == ZERO) // Return if the two lines are parallel
// The two lines are parallel: if (CGAL::sign (c[1]) == ZERO) return 0;
return (0);
xs[0] = nt_traits.convert(-c[0]) / nt_traits.convert(c[1]); xs[0] = nt_traits.convert(-c[0]) / nt_traits.convert(c[1]);
return (1); return 1;
} }
// We substitute this expression into the equation of the second // We substitute this expression into the equation of the second
// conic, and get the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0: // conic, and get the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0:
c[2] = u1*u1*s2 - u1*v1*t2 + v1*v1*r2; c[2] = u1*u1*s2 - u1*v1*t2 + v1*v1*r2;
c[1] = _two*u1*w1*s2 - u1*v1*v2 - v1*w1*t2 + v1*v1*u2; c[1] = two*u1*w1*s2 - u1*v1*v2 - v1*w1*t2 + v1*v1*u2;
c[0] = w1*w1*s2 - v1*w1*v2 + v1*v1*w2; c[0] = w1*w1*s2 - v1*w1*v2 + v1*v1*w2;
xs_end = nt_traits.solve_quadratic_equation (c[2], c[1], c[0], xs_end = nt_traits.solve_quadratic_equation(c[2], c[1], c[0], xs);
xs);
return static_cast<int>(xs_end - xs); return static_cast<int>(xs_end - xs);
} }
// At this stage, both curves have degree 2. We obtain a qaurtic polynomial // At this stage, both curves have degree 2. We obtain a qaurtic polynomial
// whose roots are the x-coordinates of the intersection points. // whose roots are the x-coordinates of the intersection points.
if (CGAL::sign (s1) == ZERO && CGAL::sign (s2) == ZERO) if (CGAL::sign (s1) == ZERO && CGAL::sign (s2) == ZERO) {
{
// If both s1 and s2 are zero, we can write the two curves as: // If both s1 and s2 are zero, we can write the two curves as:
// C1: (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0 // C1: (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0
// C2: (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0 // C2: (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0
@ -125,41 +116,40 @@ int
degree = 3; degree = 3;
} }
else else {
{
// We can write the two curves as: // We can write the two curves as:
// C1: (s1)*y^2 + (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0 // C1: (s1)*y^2 + (t1*x + v1)*y + (r1*x^2 + u1*x + w1) = 0
// C2: (s2)*y^2 + (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0 // C2: (s2)*y^2 + (t2*x + v2)*y + (r2*x^2 + u2*x + w2) = 0
// By writing the resultant of these two polynomials we get a quartic // By writing the resultant of these two polynomials we get a quartic
// polynomial, whose coefficients are given by: // polynomial, whose coefficients are given by:
c[4] = -_two*s1*s2*r1*r2 + s1*t2*t2*r1 - s1*t2*t1*r2 + c[4] = -two*s1*s2*r1*r2 + s1*t2*t2*r1 - s1*t2*t1*r2 +
s1*s1*r2*r2 - s2*t1*r1*t2 + s2*t1*t1*r2 + s2*s2*r1*r1; s1*s1*r2*r2 - s2*t1*r1*t2 + s2*t1*t1*r2 + s2*s2*r1*r1;
c[3] = -t2*r1*v1*s2 - u2*t1*t2*s1 - v2*r1*t1*s2 - c[3] = -t2*r1*v1*s2 - u2*t1*t2*s1 - v2*r1*t1*s2 -
r2*t1*v2*s1 - _two*s1*s2*r1*u2 - t2*u1*t1*s2 + u2*t1*t1*s2 - r2*t1*v2*s1 - two*s1*s2*r1*u2 - t2*u1*t1*s2 + u2*t1*t1*s2 -
r2*v1*t2*s1 + u1*t2*t2*s1 + _two*v2*r1*t2*s1 + _two*u2*r2*s1*s1 + r2*v1*t2*s1 + u1*t2*t2*s1 + two*v2*r1*t2*s1 + two*u2*r2*s1*s1 +
_two*r2*v1*t1*s2 + _two*u1*r1*s2*s2 - _two*s1*s2*u1*r2; two*r2*v1*t1*s2 + two*u1*r1*s2*s2 - two*s1*s2*u1*r2;
c[2] = -r2*v1*v2*s1 + u2*u2*s1*s1 + _two*w2*r2*s1*s1 + c[2] = -r2*v1*v2*s1 + u2*u2*s1*s1 + two*w2*r2*s1*s1 +
_two*u2*v1*t1*s2 - u2*v1*t2*s1 + w2*t1*t1*s2 - _two*s1*s2*u1*u2 - two*u2*v1*t1*s2 - u2*v1*t2*s1 + w2*t1*t1*s2 - two*s1*s2*u1*u2 -
w2*t1*t2*s1 + v2*v2*r1*s1 + u1*u1*s2*s2 - v2*r1*v1*s2 + w2*t1*t2*s1 + v2*v2*r1*s1 + u1*u1*s2*s2 - v2*r1*v1*s2 +
_two*w1*r1*s2*s2 - u2*t1*v2*s1 - t2*u1*v1*s2 - _two*s1*s2*r1*w2 - two*w1*r1*s2*s2 - u2*t1*v2*s1 - t2*u1*v1*s2 - two*s1*s2*r1*w2 -
_two*s1*s2*w1*r2 + r2*v1*v1*s2 + w1*t2*t2*s1 - v2*u1*t1*s2 - two*s1*s2*w1*r2 + r2*v1*v1*s2 + w1*t2*t2*s1 - v2*u1*t1*s2 -
t2*w1*t1*s2 + _two*v2*u1*t2*s1; t2*w1*t1*s2 + two*v2*u1*t2*s1;
c[1] = _two*w2*u2*s1*s1 + _two*w2*v1*t1*s2 - w2*v1*t2*s1 + c[1] = two*w2*u2*s1*s1 + two*w2*v1*t1*s2 - w2*v1*t2*s1 +
_two*v2*w1*t2*s1 + _two*w1*u1*s2*s2 - v2*u1*v1*s2 - _two*s1*s2*u1*w2 - two*v2*w1*t2*s1 + two*w1*u1*s2*s2 - v2*u1*v1*s2 - two*s1*s2*u1*w2 -
v2*w1*t1*s2 + u2*v1*v1*s2 - t2*w1*v1*s2 - w2*t1*v2*s1 + v2*w1*t1*s2 + u2*v1*v1*s2 - t2*w1*v1*s2 - w2*t1*v2*s1 +
v2*v2*u1*s1 - u2*v1*v2*s1 - _two*s1*s2*w1*u2; v2*v2*u1*s1 - u2*v1*v2*s1 - two*s1*s2*w1*u2;
c[0] = s2*v1*v1*w2 - s1*v2*v1*w2 - s2*v1*w1*v2 + s2*s2*w1*w1 - c[0] = s2*v1*v1*w2 - s1*v2*v1*w2 - s2*v1*w1*v2 + s2*s2*w1*w1 -
_two*s1*s2*w1*w2 + s1*w1*v2*v2 + s1*s1*w2*w2; two*s1*s2*w1*w2 + s1*w1*v2*v2 + s1*s1*w2*w2;
degree = 4; degree = 4;
} }
// Compute the roots of the resultant polynomial. // Compute the roots of the resultant polynomial.
typename Nt_traits::Polynomial poly = typename Nt_traits::Polynomial poly =
nt_traits.construct_polynomial (c, degree); nt_traits.construct_polynomial (c, degree);
xs_end = nt_traits.compute_polynomial_roots (poly, xs_end = nt_traits.compute_polynomial_roots (poly,
@ -167,8 +157,7 @@ int
return static_cast<int>(xs_end - xs); return static_cast<int>(xs_end - xs);
} }
/*! /*! Compute the roots of the resultants of the two bivariate polynomials:
* Compute the roots of the resultants of the two bivariate polynomials:
* C1: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0 * C1: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0
* C2: A*x + B*y + C = 0 * C2: A*x + B*y + C = 0
* \param deg1 The degree of the first curve. * \param deg1 The degree of the first curve.
@ -177,55 +166,50 @@ int
* \pre xs must be a vector of size 4. * \pre xs must be a vector of size 4.
* \return The number of distinct roots found. * \return The number of distinct roots found.
*/ */
template <class Nt_traits> template <typename Nt_traits>
int int compute_resultant_roots(Nt_traits& nt_traits,
_compute_resultant_roots (Nt_traits& nt_traits, const typename Nt_traits::Algebraic& r,
const typename Nt_traits::Algebraic& r, const typename Nt_traits::Algebraic& s,
const typename Nt_traits::Algebraic& s, const typename Nt_traits::Algebraic& t,
const typename Nt_traits::Algebraic& t, const typename Nt_traits::Algebraic& u,
const typename Nt_traits::Algebraic& u, const typename Nt_traits::Algebraic& v,
const typename Nt_traits::Algebraic& v, const typename Nt_traits::Algebraic& w,
const typename Nt_traits::Algebraic& w, const int& deg1,
const int& deg1, const typename Nt_traits::Algebraic& A,
const typename Nt_traits::Algebraic& A, const typename Nt_traits::Algebraic& B,
const typename Nt_traits::Algebraic& B, const typename Nt_traits::Algebraic& C,
const typename Nt_traits::Algebraic& C, typename Nt_traits::Algebraic* xs)
typename Nt_traits::Algebraic *xs)
{ {
if (deg1 == 1) if (deg1 == 1) {
{
// We should actually compute the intersection of two line: // We should actually compute the intersection of two line:
// (u*x + v*y + w = 0) and (A*x + B*y + C = 0): // (u*x + v*y + w = 0) and (A*x + B*y + C = 0):
const typename Nt_traits::Algebraic denom = A*v - B*u; const typename Nt_traits::Algebraic denom = A*v - B*u;
if (CGAL::sign (denom) == CGAL::ZERO) // Return if the two lines are parallel and do not intersect.
// The two lines are parallel and do not intersect. if (CGAL::sign(denom) == CGAL::ZERO) return 0;
return (0);
xs[0] = (B*w - C*v) / denom; xs[0] = (B*w - C*v) / denom;
return (1); return 1;
} }
if (CGAL::sign (B) == CGAL::ZERO) if (CGAL::sign(B) == CGAL::ZERO) {
{
// The first line is A*x + C = 0, therefore: // The first line is A*x + C = 0, therefore:
xs[0] = -C / A; xs[0] = -C / A;
return (1); return 1;
} }
// We can write the first curve as: y = -(A*x + C) / B. // We can write the first curve as: y = -(A*x + C) / B.
// We substitute this expression into the equation of the conic, and get // We substitute this expression into the equation of the conic, and get
// the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0: // the quadratic equation c[2]*x^2 + c[1]*x + c[0] = 0:
const typename Nt_traits::Algebraic _two = 2; const typename Nt_traits::Algebraic two = 2;
typename Nt_traits::Algebraic c[3]; typename Nt_traits::Algebraic c[3];
typename Nt_traits::Algebraic *xs_end; typename Nt_traits::Algebraic* xs_end;
c[2] = A*A*s - A*B*t + B*B*r; c[2] = A*A*s - A*B*t + B*B*r;
c[1] = _two*A*C*s - A*B*v - B*C*t + B*B*u; c[1] = two*A*C*s - A*B*v - B*C*t + B*B*u;
c[0] = C*C*s - B*C*v + B*B*w; c[0] = C*C*s - B*C*v + B*B*w;
xs_end = nt_traits.solve_quadratic_equation (c[2], c[1], c[0], xs_end = nt_traits.solve_quadratic_equation(c[2], c[1], c[0], xs);
xs);
return static_cast<int>(xs_end - xs); return static_cast<int>(xs_end - xs);
} }