WIP [scip ci]

This commit is contained in:
Andreas Fabri 2023-02-07 21:18:31 +00:00
parent fb71e3937e
commit 2184d22c52
10 changed files with 79 additions and 31 deletions

View File

@ -28,7 +28,7 @@
#include <vector> #include <vector>
#if CGAL_USE_CORE #if CGAL_USE_CORE
namespace CORE { class BigInt; } #include <CGAL/CORE/BigInt.h>
#endif #endif
namespace CGAL { namespace CGAL {

View File

@ -1389,9 +1389,9 @@ private:
// //
Nt_traits nt_traits; Nt_traits nt_traits;
const int or_fact = (_orient == CLOCKWISE) ? -1 : 1; const int or_fact = (_orient == CLOCKWISE) ? -1 : 1;
const Algebraic r = nt_traits.convert (or_fact * _r); const Algebraic r = nt_traits.convert (Integer(or_fact * _r));
const Algebraic s = nt_traits.convert (or_fact * _s); const Algebraic s = nt_traits.convert (Integer(or_fact * _s));
const Algebraic t = nt_traits.convert (or_fact * _t); const Algebraic t = nt_traits.convert (Integer(or_fact * _t));
const Algebraic cos_2phi = (r - s) / nt_traits.sqrt((r-s)*(r-s) + t*t); const Algebraic cos_2phi = (r - s) / nt_traits.sqrt((r-s)*(r-s) + t*t);
const Algebraic _zero = 0; const Algebraic _zero = 0;
const Algebraic _one = 1; const Algebraic _one = 1;
@ -1441,8 +1441,8 @@ private:
// 4*r*s - t^2 4*r*s - t^2 // 4*r*s - t^2 4*r*s - t^2
// //
// The denominator (4*r*s - t^2) must be negative for hyperbolas. // The denominator (4*r*s - t^2) must be negative for hyperbolas.
const Algebraic u = nt_traits.convert (or_fact * _u); const Algebraic u = nt_traits.convert (Integer(or_fact * _u));
const Algebraic v = nt_traits.convert (or_fact * _v); const Algebraic v = nt_traits.convert (Integer(or_fact * _v));
const Algebraic det = 4*r*s - t*t; const Algebraic det = 4*r*s - t*t;
Algebraic x0, y0; Algebraic x0, y0;
@ -1650,9 +1650,9 @@ protected:
int n_xs; int n_xs;
Nt_traits nt_traits; Nt_traits nt_traits;
xs_end = nt_traits.solve_quadratic_equation (_t*_t - _four*_r*_s, xs_end = nt_traits.solve_quadratic_equation (Integer(_t*_t - _four*_r*_s),
_two*_t*_v - _four*_s*_u, Integer(_two*_t*_v - _four*_s*_u),
_v*_v - _four*_s*_w, Integer(_v*_v - _four*_s*_w),
xs); xs);
n_xs = static_cast<int>(xs_end - xs); n_xs = static_cast<int>(xs_end - xs);
@ -1664,15 +1664,15 @@ protected:
if (CGAL::sign (_t) == ZERO) if (CGAL::sign (_t) == ZERO)
{ {
// The two vertical tangency points have the same y coordinate: // The two vertical tangency points have the same y coordinate:
ys[0] = nt_traits.convert (-_v) /nt_traits.convert (_two*_s); ys[0] = nt_traits.convert (Integer(-_v)) /nt_traits.convert (Integer(_two*_s));
n_ys = 1; n_ys = 1;
} }
else else
{ {
ys_end = ys_end =
nt_traits.solve_quadratic_equation (_four*_r*_s*_s - _s*_t*_t, nt_traits.solve_quadratic_equation (Integer(_four*_r*_s*_s - _s*_t*_t),
_four*_r*_s*_v - _two*_s*_t*_u, Integer(_four*_r*_s*_v - _two*_s*_t*_u),
_r*_v*_v - _t*_u*_v + _t*_t*_w, Integer(_r*_v*_v - _t*_u*_v + _t*_t*_w),
ys); ys);
n_ys = static_cast<int>(ys_end - ys); n_ys = static_cast<int>(ys_end - ys);
} }
@ -1692,7 +1692,7 @@ protected:
{ {
for (j = 0; j < n_ys; j++) for (j = 0; j < n_ys; j++)
{ {
if (CGAL::compare (nt_traits.convert(_two*_s) * ys[j], if (CGAL::compare (nt_traits.convert(Integer(_two*_s)) * ys[j],
-(nt_traits.convert(_t) * xs[i] + -(nt_traits.convert(_t) * xs[i] +
nt_traits.convert(_v))) == EQUAL) nt_traits.convert(_v))) == EQUAL)
{ {
@ -1735,9 +1735,9 @@ protected:
Algebraic *ys_end; Algebraic *ys_end;
Nt_traits nt_traits; Nt_traits nt_traits;
ys_end = nt_traits.solve_quadratic_equation (_t*_t - _four*_r*_s, ys_end = nt_traits.solve_quadratic_equation (Integer(_t*_t - _four*_r*_s),
_two*_t*_u - _four*_r*_v, Integer(_two*_t*_u - _four*_r*_v),
_u*_u - _four*_r*_w, Integer(_u*_u - _four*_r*_w),
ys); ys);
n = static_cast<int>(ys_end - ys); n = static_cast<int>(ys_end - ys);

View File

@ -72,13 +72,14 @@ int
unsigned int degree = 4; unsigned int degree = 4;
typename Nt_traits::Algebraic *xs_end; typename Nt_traits::Algebraic *xs_end;
typedef typename Nt_traits::Integer Integer;
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(Integer(-w1)) / nt_traits.convert(u1);
return (1); return (1);
} }
@ -94,7 +95,7 @@ int
// The two lines are parallel: // The two lines are parallel:
return (0); return (0);
xs[0] = nt_traits.convert(-c[0]) / nt_traits.convert(c[1]); xs[0] = nt_traits.convert(Integer(-c[0])) / nt_traits.convert(c[1]);
return (1); return (1);
} }

View File

@ -225,11 +225,12 @@ public:
if (CGAL::compare( if (CGAL::compare(
CGAL::width(y_bfi), CGAL::width(y_bfi),
CGAL::lower(CGAL::abs(y_bfi)) * eps) CGAL::lower(CGAL::abs(y_bfi)) * eps)
== SMALLER) == SMALLER){
return std::make_pair( return std::make_pair(
Bound(CGAL::lower(y_bfi)), Bound(CGAL::lower(y_bfi)),
Bound(CGAL::upper(y_bfi))); Bound(CGAL::upper(y_bfi)));
} }
}
else precision*=2; else precision*=2;
} }
} }
@ -287,10 +288,12 @@ private:
if (CGAL::zero_in(y_denom_bfi) == false) if (CGAL::zero_in(y_denom_bfi) == false)
{ {
BFI y_bfi(y_numer_bfi/y_denom_bfi); BFI y_bfi(y_numer_bfi/y_denom_bfi);
if (CGAL::width(y_bfi) < eps ) if (CGAL::width(y_bfi) < eps ){
assert(false) ; // AF fix
return std::make_pair( return std::make_pair(
Bound(CGAL::lower(y_bfi)), Bound(CGAL::lower(y_bfi)),
Bound(CGAL::upper(y_bfi))); Bound(CGAL::upper(y_bfi)));
}
} }
else precision*=2; else precision*=2;

View File

@ -246,8 +246,8 @@ public:
if (sign_disc == ZERO) if (sign_disc == ZERO)
{ {
// We have one real root with mutliplicity 2. // We have one real root with multiplicity 2.
*oi = -Algebraic (b) / Algebraic (2*a); *oi = -Algebraic (b) / Algebraic (NT(2*a));
++oi; ++oi;
} }
else if (sign_disc == POSITIVE) else if (sign_disc == POSITIVE)
@ -255,7 +255,7 @@ public:
// We have two distinct real roots. We return them in ascending order. // We have two distinct real roots. We return them in ascending order.
const Algebraic sqrt_disc = CGAL::sqrt (Algebraic (disc)); const Algebraic sqrt_disc = CGAL::sqrt (Algebraic (disc));
const Algebraic alg_b = b; const Algebraic alg_b = b;
const Algebraic alg_2a = 2*a; const Algebraic alg_2a = NT(2*a);
if (sign_a == POSITIVE) if (sign_a == POSITIVE)
{ {

View File

@ -614,7 +614,10 @@ inline BigFloat gcd(const BigFloat& a, const BigFloat& b) {
//mpz_tdiv_qr(q.get_mp(), r.get_mp(), a.get_mp(), b.get_mp()); //mpz_tdiv_qr(q.get_mp(), r.get_mp(), a.get_mp(), b.get_mp());
//}// //}//
/* AF /*
// AF: As BigRat is just a boost::mp type we cannot have this
// constructor
// constructor BigRat from BigFloat // constructor BigRat from BigFloat
inline BigRat::BigRat(const BigFloat& f) : RCBigRat(new BigRatRep()){ inline BigRat::BigRat(const BigFloat& f) : RCBigRat(new BigRatRep()){
*this = f.BigRatValue(); *this = f.BigRatValue();

View File

@ -151,6 +151,9 @@ inline long div_exact(long x, long y) {
return x/y; // precondition: isDivisible(x,y) return x/y; // precondition: isDivisible(x,y)
} }
inline BigInt gcd(const BigInt& a, const BigInt& b){
return boost::multiprecision::gcd(a,b);
}
/// ceilLg -- ceiling of log_2(a) where a=BigInt, int or long /// ceilLg -- ceiling of log_2(a) where a=BigInt, int or long
/** Convention: a=0, ceilLg(a) returns -1. /** Convention: a=0, ceilLg(a) returns -1.

View File

@ -44,6 +44,43 @@ namespace CORE {
{ {
return boost::multiprecision::denominator(q); return boost::multiprecision::denominator(q);
} }
// Chee (3/19/2004):
// The following definitions of div_exact(x,y) and gcd(x,y)
// ensures that in Polynomial<NT>
/// divisible(x,y) = "x | y"
inline BigRat div_exact(const BigRat& x, const BigRat& y) {
BigRat z = x / y;
return z;
}
inline BigRat gcd(const BigRat& x , const BigRat& y)
{
// return BigRat(1); // Remark: we may want replace this by
// the definition of gcd of a quotient field
// of a UFD [Yap's book, Chap.3]
//Here is one possible definition: gcd of x and y is just the
//gcd of the numerators of x and y divided by the gcd of the
//denominators of x and y.
BigInt n = gcd(numerator(x), numerator(y));
BigInt d = gcd(denominator(x), denominator(y));
return BigRat(n,d);
}
// Chee: 8/8/2004: need isDivisible to compile Polynomial<BigRat>
// A trivial implementation is to return true always. But this
// caused tPolyRat to fail.
// So we follow the definition of
// Expr::isDivisible(e1, e2) which checks if e1/e2 is an integer.
inline bool isInteger(const BigRat& x) {
return denominator(x) == 1; // AF: does that need canonicalize?
}
inline bool isDivisible(const BigRat& x, const BigRat& y) {
BigRat r;
r = x/y;
return isInteger(r);
}
/// BigIntValue /// BigIntValue
inline BigInt BigIntValue(const BigRat& br) inline BigInt BigIntValue(const BigRat& br)
{ {

View File

@ -48,6 +48,7 @@
#define CORE_POLY_H #define CORE_POLY_H
#include <CGAL/CORE/BigFloat.h> #include <CGAL/CORE/BigFloat.h>
#include <CGAL/CORE/BigRat.h>
#include <CGAL/CORE/Promote.h> #include <CGAL/CORE/Promote.h>
#include <vector> #include <vector>
#include <CGAL/assertions.h> #include <CGAL/assertions.h>

View File

@ -892,10 +892,10 @@ BigFloat Polynomial<NT>::CauchyUpperBound() const {
NT mx = 0; NT mx = 0;
int deg = getTrueDegree(); int deg = getTrueDegree();
for (int i = 0; i < deg; ++i) { for (int i = 0; i < deg; ++i) {
mx = core_max(mx, abs(coeff[i])); mx = core_max(mx, NT(abs(coeff[i])));
} }
Expr e = mx; Expr e = mx;
e /= Expr(abs(coeff[deg])); e /= Expr(NT(abs(coeff[deg])));
e.approx(CORE_INFTY, 2); e.approx(CORE_INFTY, 2);
// get an absolute approximate value with error < 1/4 // get an absolute approximate value with error < 1/4
return (e.BigFloatValue().makeExact() + 2); return (e.BigFloatValue().makeExact() + 2);
@ -975,9 +975,9 @@ BigFloat Polynomial<NT>::CauchyLowerBound() const {
NT mx = 0; NT mx = 0;
int deg = getTrueDegree(); int deg = getTrueDegree();
for (int i = 1; i <= deg; ++i) { for (int i = 1; i <= deg; ++i) {
mx = core_max(mx, abs(coeff[i])); mx = core_max(mx, NT(abs(coeff[i])));
} }
Expr e = Expr(abs(coeff[0]))/ Expr(abs(coeff[0]) + mx); Expr e = Expr(NT(abs(coeff[0])))/ Expr(NT(abs(coeff[0])) + mx);
e.approx(2, CORE_INFTY); e.approx(2, CORE_INFTY);
// get an relative approximate value with error < 1/4 // get an relative approximate value with error < 1/4
return (e.BigFloatValue().makeExact().div2()); return (e.BigFloatValue().makeExact().div2());
@ -1018,8 +1018,8 @@ BigFloat Polynomial<NT>::height() const {
int deg = getTrueDegree(); int deg = getTrueDegree();
NT ht = 0; NT ht = 0;
for (int i = 0; i< deg; i++) for (int i = 0; i< deg; i++)
if (ht < abs(coeff[i])) if (ht < NT(abs(coeff[i])))
ht = abs(coeff[i]); ht = NT(abs(coeff[i]));
return BigFloat(ht); return BigFloat(ht);
} }