songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Polynomial/test/Polynomial/Polynomial.cpp

1198 lines
38 KiB
C++
Raw Blame History

#include <CGAL/basic.h>
#include <CGAL/Test/_test_algebraic_structure.h>
#include <CGAL/Test/_test_real_embeddable.h>
#include <iostream>
#include <cassert>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Arithmetic_kernel.h>
#include <CGAL/Polynomial.h>
#include <CGAL/Sqrt_extension.h>
#include <CGAL/Interval_nt.h>
#include <CGAL/Coercion_traits.h>
#define CGAL_DEFINE_TYPES_FROM_AK(AK) \
typedef typename AK::Integer Integer; \
typedef typename AK::Rational Rational; \
typedef typename AK::Field_with_sqrt Field_with_sqrt; \
typedef CGAL::Polynomial<Integer> Poly_int1; \
typedef CGAL::Polynomial<Poly_int1> Poly_int2; \
typedef CGAL::Polynomial<Poly_int2> Poly_int3; \
typedef CGAL::Polynomial<Rational> Poly_rat1; \
typedef CGAL::Polynomial<Poly_rat1> Poly_rat2; \
typedef CGAL::Polynomial<Poly_rat2> Poly_rat3;
// TODO: copied from number_type_utils.h
CGAL_BEGIN_NAMESPACE
template <class NT , class RT>
inline
void convert_to(const NT& x, RT& r){
typedef CGAL::Coercion_traits<NT,RT> CT;
typedef typename CT::Coercion_type RET;
BOOST_STATIC_ASSERT((::boost::is_same<RET,RT>::value));
r = typename CT::Cast()(x);
}
CGAL_END_NAMESPACE
template <class NT>
void datastru() {
typedef CGAL::Polynomial<NT> POLY;
int i;
// uninitialized construction
POLY p;
assert( p.degree() == 0 );
// construction from iterator range, assignment
NT array[] = { NT(3), NT(1), NT(4), NT(1), NT(5), NT(9) };
NT* const array_end = array + sizeof array / sizeof *array;
POLY q(array, array_end);
p = q;
assert( array_end - array == p.degree() + 1 );
assert( std::equal(p.begin(), p.end(), array) );
// immediate construction, reduce(), copy construction, (in)equality
POLY r(NT(3), NT(1), NT(4), NT(1), NT(5), NT(9), NT(0), NT(0));
assert( q.degree() == r.degree() );
assert( q == r );
assert( r == r );
q = POLY(NT(2), NT(7), NT(1), NT(8), NT(1));
assert( q != r );
POLY s(p);
assert( r == s );
// selection of coefficients
for (i = 0; i <= p.degree(); i++) {
assert( p[i] == array[i] );
}
assert( i == array_end - array );
assert( p.lcoeff() == array_end[-1] );
// reversal()
assert( reversal(q) == POLY(NT(1), NT(8), NT(1), NT(7), NT(2)) );
// divide_by_x()
r.divide_by_x();
assert( r.degree() == 4 );
assert( r == POLY(NT(1), NT(4), NT(1), NT(5), NT(9)) );
// zero test
assert( !p.is_zero() );
q = POLY(NT(0));
assert( q.is_zero() );
}
template <class NT> void signs(::CGAL::Tag_false);
template <class NT> void signs(::CGAL::Tag_true);
template <class NT> void division(CGAL::Integral_domain_without_division_tag);
template <class NT> void division(CGAL::Integral_domain_tag);
template <class NT>
void arithmetic() {
typedef CGAL::Polynomial<NT> POLY;
// evaluation
POLY p(NT(6), NT(-5), NT(1));
assert( p.evaluate(NT(1)) == NT(2) );
assert( p.evaluate(NT(2)) == NT(0) );
assert( p.evaluate(NT(3)) == NT(0) );
// TODO: No evaluate_absolute available for all NTs.
//assert( p.evaluate_absolute(NT(1)) == NT(12) );
// diff
POLY q(NT(3), NT(2), NT(-2));
POLY z(NT(0));
assert( diff(z) == z );
assert( p.degree() == 2 );
assert( diff(p) == POLY(NT(-5), NT(2)) );
assert( p.degree() == 2 );
assert( q.degree() == 2 );
q.diff();
assert( q == POLY(NT(2), NT(-4)) );
// scale and translate
p.scale_up(NT(4));
assert( p == POLY(NT(6), NT(-20), NT(16)) );
assert( scale_down(p, NT(2)) == POLY(NT(24), NT(-40), NT(16)) );
assert( p == POLY(NT(6), NT(-20), NT(16)) );
assert( translate_by_one(p) == translate(p, NT(1)) );
p = POLY(NT(3), NT(2), NT(1));
p.translate(NT(-2));
assert( p == POLY(NT(3), NT(-2), NT(1)) );
// addition
POLY r;
p = POLY(NT(3), NT(2), NT(1), NT(1));
q = POLY(NT(1), NT(2), NT(-1), NT(-1));
r = p + q;
assert( r.degree() == 1 );
assert( r[0] == NT(4) );
assert( r[1] == NT(4) );
q += p;
assert( q == r );
q += NT(1);
assert( q == POLY(NT(5), NT(4)) );
assert( q == r + NT(1) );
assert( q == NT(1) + r );
// subtraction
p = POLY(NT(3), NT(2), NT(1), NT(1));
q = POLY(NT(1), NT(2), NT(-1), NT(1));
r = p - q;
assert( r.degree() == 2 );
assert( r == POLY(NT(2), NT(0), NT(2)) );
r -= -q;
assert( r == p );
r -= NT(2);
assert( r == p - NT(2) );
assert( -r == NT(2) - p );
// multiplication
p = POLY(NT(1), NT(2), NT(3));
q = POLY(NT(-2), NT(5), NT(2));
r = p * q;
assert( r == POLY(NT(-2), NT(1), NT(6), NT(19), NT(6)) );
q *= p;
assert( q == r );
q *= NT(2);
assert( q == POLY(NT(-4), NT(2), NT(12), NT(38), NT(12)) );
assert( q == r * NT(2) );
assert( q == NT(2) * r );
// sign etc., if applicable
// TODO: Replaced Is_real_comparable with RET::Is_real_embeddable, OK?
typedef typename CGAL::Real_embeddable_traits<NT>::Is_real_embeddable
Embeddable;
signs<NT>(Embeddable());
// division, if applicable
typedef typename CGAL::Algebraic_structure_traits<NT>::Algebraic_category
Algebra_type;
division<NT>(Algebra_type());
}
template <class NT> void signs(::CGAL::Tag_false) {
// nothing to be done
}
template <class POLY>
void test_greater(POLY a, POLY b) {
// check .compare()
assert( a.compare(b) == CGAL::LARGER );
assert( a.compare(a) == CGAL::EQUAL );
assert( b.compare(a) == CGAL::SMALLER );
// check that all comparison operators reflect a > b
assert( a > b );
assert( !(b > a) );
assert( !(a > a) );
assert( b < a );
assert( !(a < b) );
assert( !(a < a) );
assert( a >= b );
assert( !(b >= a) );
assert( a >= a );
assert( b <= a );
assert( !(a <= b) );
assert( a <= a );
}
template <class NT>
void signs(::CGAL::Tag_true) {
typedef CGAL::Polynomial<NT> POLY;
POLY p(NT(6), NT(-5), NT(1));
POLY pp(NT(6), NT(-5), NT(1));
POLY q(NT(3), NT(2), NT(-2));
POLY z(NT(0));
// sign, abs
assert( z.sign() == CGAL::ZERO );
assert( p.sign() == CGAL::POSITIVE );
assert( p.abs() == p );
assert( q.sign() == CGAL::NEGATIVE );
assert( q.abs() == POLY(NT(-3), NT(-2), NT(2)) );
// comparison operators
test_greater(p + NT(1), p);
test_greater(POLY(NT(1), NT(2), NT(3), NT(4)), p);
// sign_at
q=POLY(NT(3),NT(2),NT(-1));
assert(q.sign_at(-2)==CGAL::NEGATIVE);
assert(q.sign_at(-1)==CGAL::ZERO);
assert(q.sign_at( 0)==CGAL::POSITIVE);
assert(q.sign_at( 3)==CGAL::ZERO);
assert(q.sign_at( 4)==CGAL::NEGATIVE);
assert(z.sign_at( 3)==CGAL::ZERO);
}
template <class NT>
void division(CGAL::Integral_domain_without_division_tag) {
// nothing to be done
}
template <class NT>
void division(CGAL::Integral_domain_tag) {
typedef CGAL::Polynomial<NT> POLY;
// integral division (remainder zero)
POLY p(NT(2), NT(4), NT(6)), q(NT(-2), NT(5), NT(2));
POLY r(NT(-4), NT(2), NT(12), NT(38), NT(12)), s;
assert( q == r / p );
r /= q;
assert( r == p );
r /= NT(2);
assert( r == POLY(NT(1), NT(2), NT(3)) );
assert( NT(6) / POLY(NT(3)) == POLY(NT(2)) );
// division with remainder
p = POLY(NT(2), NT(-3), NT(1)); // (x-1)(x-2)
q = POLY(NT(24), NT(-14), NT(2)); // 2(x-3)(x-4)
r = POLY(NT(-5), NT(1)); // (x-5)
s = p*q + r;
assert(!CGAL::divides(q,s));
POLY Q, R;
POLY::euclidean_division(s, p, Q, R); // works, since p.lcoeff() == 1
assert( Q == q );
assert( R == r );
NT D;
POLY::pseudo_division(s, q, Q, R, D);
assert( Q == D*p );
assert( R == D*r );
// codecover
p = POLY(NT(2),NT(3));
s = POLY(NT(2),NT(3),NT(4));
POLY::euclidean_division(p, s, Q, R); // works, since p.lcoeff() == 1
assert( Q == POLY(NT(0)) );
assert( R == p );
POLY::pseudo_division(p, s, Q, R, D); // works, since p.lcoeff() == 1
assert( Q == POLY(NT(0)) );
assert( R == D*p );
}
template <class NT>
void io() {
typedef CGAL::Polynomial<NT> POLY;
{
// successful re-reading of output
POLY p(NT(-3), NT(5), NT(0), NT(0), NT(-7), NT(9)), q;
std::ostringstream os;
os << p;
std::istringstream is(os.str());
is >> q;
assert( p == q );
}{
std::ostringstream os;
CGAL::set_pretty_mode(os);
os << oformat(POLY(NT(3)));
//std::cout <<os.str()<<std::endl;
assert( os.str() == "3" );
}{
std::ostringstream os;
CGAL::set_pretty_mode(os);
os << oformat(POLY(NT(-3)));
assert( os.str() == "(-3)" );
}{
std::ostringstream os;
CGAL::set_pretty_mode(os);
os << oformat(POLY(NT(-3)),CGAL::Parens_as_product_tag());
assert( os.str() == "(-3)" );
}{
std::ostringstream os;
CGAL::set_pretty_mode(os);
os << oformat(POLY(NT(-3),NT(4)));
if( CGAL::Polynomial_traits_d<POLY>::d == 1)
assert( os.str() == "4*x + (-3)" );
else
assert( os.str() == "4*y + (-3)" );
}{
std::ostringstream os;
CGAL::set_pretty_mode(os);
os << oformat(POLY(NT(-3),NT(4)), CGAL::Parens_as_product_tag());
if( CGAL::Polynomial_traits_d<POLY>::d == 1)
assert( os.str() == "(4*x + (-3))" );
else
assert( os.str() == "(4*y + (-3))" );
}
}
template <class NT>
void canon(CGAL::Integral_domain_without_division_tag) {
// dummy for cases where canonicalize_polynomial() does not apply
}
template <class NT>
void canon(CGAL::Unique_factorization_domain_tag) {
typedef CGAL::Polynomial<NT> UNPOLY;
typedef CGAL::Polynomial<UNPOLY> BIPOLY;
BIPOLY p(0), c;
assert(canonicalize_polynomial(p) == BIPOLY(0));
p = BIPOLY(UNPOLY(NT(15)), UNPOLY(NT(-12), NT( 9)), UNPOLY(NT(0), NT(-6)));
c = BIPOLY(UNPOLY(NT(-5)), UNPOLY(NT( 4), NT(-3)), UNPOLY(NT(0), NT( 2)));
assert(canonicalize_polynomial(p) == c);
}
template <class NT>
void canon(CGAL::Field_tag) {
typedef CGAL::Polynomial<NT> UNPOLY;
typedef CGAL::Polynomial<UNPOLY> BIPOLY;
BIPOLY p(0), c;
assert(canonicalize_polynomial(p) == BIPOLY(0));
p = BIPOLY(UNPOLY(NT(15)), UNPOLY(NT(-12), NT( 9)), UNPOLY(NT(0), NT(-6)));
c = BIPOLY(UNPOLY(NT(-5)/NT(2)), UNPOLY(NT(2), NT(-3)/NT(2)),
UNPOLY(NT(0), NT(1)));
assert(canonicalize_polynomial(p) == c);
}
#ifdef NiX_POLY_USE_NT_TESTS
template <class NT>
void nt_tests(::CGAL::Tag_false) {
typedef CGAL::Polynomial<NT> POLY;
if ((volatile bool)0) { // test compilation only
CGAL::test_number_type<POLY,
typename CGAL::Algebraic_structure_traits<POLY>::Algebric_structure_tag>();
}
}
template <class NT>
void nt_tests(::CGAL::Tag_true) {
typedef CGAL::Polynomial<NT> POLY;
CGAL::test_real_comparable<POLY>();
nt_tests<NT>(::CGAL::Tag_false());
}
#endif // NiX_POLY_USE_NT_TESTS
template <class NT>
void unigcdres(CGAL::Field_tag) {
typedef CGAL::Polynomial<NT> POLY;
// test univariate polynomial gcd and resultant over a field
POLY f(NT(2), NT(7), NT(1), NT(8), NT(1), NT(8));
POLY g(NT(3), NT(1), NT(4), NT(1), NT(5), NT(9));
POLY h(NT(3), NT(4), NT(7), NT(7));
f /= NT(3);
g /= NT(5);
h /= h.unit_part();
POLY d;
d = gcd(f, g);
assert( d == POLY(NT(1)) );
assert( prs_resultant(f, g) == NT(230664271L)/NT(759375L) ); // Maple
POLY fh(f*h), gh(g*h);
d = gcd(fh, gh);
assert( d == h );
assert( prs_resultant(fh, gh) == NT(0) );
POLY a, b;
d = gcdex(fh, gh, a, b);
assert( d == h );
assert( d == a*fh + b*gh );
}
template <class NT>
void unigcdres(CGAL::Integral_domain_tag) {
typedef CGAL::Polynomial<NT> POLY;
// test univariate polynomial gcd and resultant over a non-field
POLY f(NT(2), NT(7), NT(1), NT(8), NT(1), NT(8));
POLY g(NT(3), NT(1), NT(4), NT(1), NT(5), NT(9));
POLY h(NT(3), NT(4), NT(7), NT(7));
POLY d;
NT c(42);
d = gcd((-c)*f, c*g);
assert( d == POLY(c) );
assert( prs_resultant(f, g) == NT(230664271L) ); // as computed by Maple
POLY fh(f*h), gh(g*h);
d = gcd((-c)*fh, c*gh);
assert( d == c*h );
assert( prs_resultant(fh, gh) == NT(0) );
POLY a, b; NT v;
d = pseudo_gcdex((-c)*fh, c*gh, a, b, v);
assert( d == c*h );
assert( v*d == (-c)*a*fh + c*b*gh );
// Michael Kerber's example for the hgdelta_update() bug:
// These polynomials cretate a situation where h does not divide g,
// but h^(delta-1) divides g^delta (as predicted by subresultant theory).
// TODO: Uncomment following code
/*CGAL::Creator_1<int, NT> int2nt;
static const int cp[] = { 196, 0, -140, 0, 69, 0, -8, 0, -4 },
cq[] = { -15156, 0, -58572, 0, -82309, 0, -27032, 0,
20190, 0, 1696, 0, -6697, 0, 1644, 0, 1664, 0, 256 };
static const size_t np = sizeof cp / sizeof *cp,
nq = sizeof cq / sizeof *cq;
POLY p(LiS::mapper(cp, int2nt), LiS::mapper(cp+np, int2nt)),
q(LiS::mapper(cq, int2nt), LiS::mapper(cq+nq, int2nt));
d = gcd(p,q);
assert ( d == POLY(1) );*/
}
template <class NT>
void bigcdres(CGAL::Field_tag) {
typedef CGAL::Polynomial<NT> POLY1;
typedef CGAL::Polynomial<POLY1> POLY2;
// test bivariate polynomial gcd and resultant over a field
POLY2 f = POLY2(POLY1(NT(2), NT(7), NT(1)),
POLY1(NT(8), NT(1)), POLY1(NT(8)));
POLY2 g = POLY2(POLY1(NT(3), NT(1), NT(4)),
POLY1(NT(1), NT(5)), POLY1(NT(9)));
POLY2 h = POLY2(POLY1(NT(3), NT(4), NT(1)),
POLY1(NT(9), NT(7)), POLY1(NT(7)));
f /= NT(3);
g /= NT(5);
h /= h.unit_part();
POLY2 d;
d = gcd(f, g);
assert( d == POLY2(1) );
POLY1 r(NT(1444), NT(-1726), NT(3295), NT(-2501), NT(560));
r /= NT(225);
assert( prs_resultant(f, g) == r ); // says Maple
POLY2 fh(f*h), gh(g*h);
d = gcd(fh, gh);
assert( d == h );
assert( prs_resultant(fh, gh) == POLY1(0) );
}
template <class NT>
void bigcdres(CGAL::Integral_domain_tag) {
typedef CGAL::Polynomial<NT> POLY1;
typedef CGAL::Polynomial<POLY1> POLY2;
// test bivariate polynomial gcd and resultant over a non-field
POLY2 f = POLY2(POLY1(NT(2), NT(7), NT(1)),
POLY1(NT(8), NT(1)), POLY1(NT(8)));
POLY2 g = POLY2(POLY1(NT(3), NT(1), NT(4)),
POLY1(NT(1), NT(5)), POLY1(NT(9)));
POLY2 h = POLY2(POLY1(NT(3), NT(4), NT(1)),
POLY1(NT(9), NT(7)), POLY1(NT(7)));
POLY2 c(42), d;
d = gcd(-c*f, c*g);
assert( d == c );
POLY1 r(NT(1444), NT(-1726), NT(3295), NT(-2501), NT(560));
assert( prs_resultant(f, g) == r ); // says Maple
POLY2 fh(f*h), gh(g*h);
d = gcd(-c*fh, c*gh);
assert( d == c*h );
assert( prs_resultant(fh, gh) == POLY1(0) );
}
template <class POLY>
void test_sqff_utcf_(const POLY& poly, int n){
typedef CGAL::Polynomial_traits_d<POLY> PT;
std::vector<POLY> fac;
std::vector<int> mul;
std::back_insert_iterator<std::vector<POLY> > fac_bi(fac);
std::back_insert_iterator<std::vector<int> > mul_bi(mul);
assert(n ==
CGAL::CGALi::square_free_factorize_utcf(poly, fac_bi, mul_bi));
assert((int) mul.size() == n);
assert((int) fac.size() == n);
if (typename PT::Total_degree()(poly) == 0 ){
assert(n == 0);
}else{
POLY check(1);
for(int i = 0; i < n ; i++){
assert(typename PT::Total_degree()(fac[i]) > 0);
//std::cout << fac[i] << std::endl;
//std::cout << typename PT::Canonicalize()(fac[i]) << std::endl;
//std::cout << std::endl;
assert(fac[i] == typename PT::Canonicalize()(fac[i]));
check *= CGAL::ipower(fac[i],mul[i]);
}
assert(check == typename PT::Canonicalize()(poly));
}
}
template <class POLY_1>
void test_sqff_utcf(const POLY_1& a1,const POLY_1& b1,const POLY_1& c1){
typedef CGAL::Polynomial<POLY_1> POLY_2;
test_sqff_utcf_(POLY_1(0), 0);
test_sqff_utcf_(POLY_1(1), 0);
test_sqff_utcf_(POLY_1(2), 0);
test_sqff_utcf_(a1, 1);
test_sqff_utcf_(a1*a1, 1);
test_sqff_utcf_(a1*a1*b1,2);
test_sqff_utcf_(a1*POLY_1(5),1);
POLY_2 a2(a1,b1,c1);
POLY_2 b2(a1,c1);
POLY_2 c2(b1,a1,b1);
test_sqff_utcf_(POLY_2(0), 0);
test_sqff_utcf_(POLY_2(1), 0);
test_sqff_utcf_(POLY_2(2), 0);
test_sqff_utcf_(a2,1);
test_sqff_utcf_(a2*a2,1);
test_sqff_utcf_(a2*a2*b2,2);
test_sqff_utcf_(a2*POLY_2(5),1);
test_sqff_utcf_(a2*a2*POLY_2(5),1);
// non regular polynomials
test_sqff_utcf_(a1*a2,2);
test_sqff_utcf_(a1*a2*a2,2);
test_sqff_utcf_(a1*a2*b2*b2,3);
test_sqff_utcf_(a1*a2*POLY_2(5),2);
test_sqff_utcf_(a1*a2*a2*POLY_2(5),2);
test_sqff_utcf_(a1*a2*a2*b2*POLY_2(5),3);
test_sqff_utcf_(a1*a1*a2,2);
test_sqff_utcf_(a1*a1*a2*a2,2);
test_sqff_utcf_(a1*b1*b1*a2*a2,3);
}
template <class AT>
void psqff(){
{
typedef typename AT::Integer NT;
typedef CGAL::Polynomial<NT> POLY_1;
// monic factors
POLY_1 a1(NT("4352435"), NT("4325245"));
POLY_1 b1(NT("123"), NT("432235"),NT("43324252"));
POLY_1 c1(NT("12324"), NT("25332235"),NT("24657252"));
test_sqff_utcf(a1,b1,c1);
}{
typedef typename AT::Rational NT;
typedef CGAL::Polynomial<NT> POLY_1;
// monic factors
POLY_1 a1(NT("4352435"), NT("4325245"));
POLY_1 b1(NT("123"), NT("432235"),NT("43324252"));
POLY_1 c1(NT("12324"), NT("25332235"),NT("24657252"));
test_sqff_utcf(a1,b1,c1);
}
{
typedef typename AT::Integer Integer; //UFD domain
typedef CGAL::Sqrt_extension<Integer,Integer> NT;
typedef CGAL::Polynomial<NT> POLY;
// square-free factorization (i.e. factorization by multiplicities)
POLY a1(NT(3,4,7), NT(2,3,7), NT(3));
POLY b1(NT(9,2,7), NT(5,11,7), NT(2));
POLY c1(NT(2,8,7), NT(6,3,7), NT(5));
test_sqff_utcf(a1,b1,c1);
}{
typedef typename AT::Integer Integer; //UFD domain
typedef typename AT::Rational Rational; //UFD domain
typedef CGAL::Sqrt_extension<Rational,Integer> NT;
typedef CGAL::Polynomial<NT> POLY;
// square-free factorization (i.e. factorization by multiplicities)
POLY a1(NT(3,4,7), NT(2,3,7), NT(3));
POLY b1(NT(9,2,7), NT(5,11,7), NT(2));
POLY c1(NT(2,8,7), NT(6,3,7), NT(5));
test_sqff_utcf(a1,b1,c1);
}
{
typedef typename AT::Integer Integer;
typedef typename AT::Rational Rational;
typedef CGAL::Sqrt_extension<Rational,Rational> NT;
typedef CGAL::Polynomial<NT> POLY;
// square-free factorization (i.e. factorization by multiplicities)
POLY a1(NT(3,4,7), NT(2,3,7), NT(3));
POLY b1(NT(9,2,7), NT(5,11,7), NT(2));
POLY c1(NT(2,8,7), NT(6,3,7), NT(5));
test_sqff_utcf(a1,b1,c1);
}
}
template <class NT>
void sqff() {
typedef CGAL::Polynomial<NT> POLY;
// square-free factorization (i.e. factorization by multiplicities)
POLY p1(NT(3), NT(4), NT(1));
POLY p3(NT(9), NT(5), NT(1));
POLY p4(NT(2), NT(6), NT(1));
POLY p = NT(3) * p1 * p3*p3*p3 * p4*p4*p4*p4;
typedef std::vector<POLY> PVEC;
typedef std::vector<int> IVEC;
PVEC fac;
IVEC mul;
std::back_insert_iterator<PVEC> fac_bi(fac);
std::back_insert_iterator<IVEC> mul_bi(mul);
unsigned n;
NT alpha;
n = CGAL::CGALi::square_free_factorize(p, fac_bi, mul_bi );
// assert(alpha == 3);
assert(n == 3);
assert(mul.size() == n);
assert(fac.size() == n);
assert(mul[0] == 1 && fac[0] == p1);
assert(mul[1] == 3 && fac[1] == p3);
assert(mul[2] == 4 && fac[2] == p4);
/*p = POLY( NT(1), NT(-2), NT(1) );
std::cerr << p << std::endl;
fac.clear();
mul.clear();
fac_bi = std::back_insert_iterator<PVEC>(fac);
mul_bi = std::back_insert_iterator<IVEC>(mul);
std::cerr << CGAL::CGALi::square_free_factorize( p, fac_bi, mul_bi ) << std::endl;
std::cerr << fac[0] << std::endl;*/
//std::cerr << fac[1] << std::endl;
typedef CGAL::Polynomial< CGAL::Polynomial< NT > > BPOLY;
typedef std::vector< BPOLY > BPVEC;
BPOLY bp( POLY( NT(1), NT(-2), NT(1) ) );
// std::cerr << bp << std::endl;
// typename CGAL::Polynomial_traits_d< BPOLY >::Make_square_free make_square_free;
// std::cerr << make_square_free( bp ) << std::endl;
// std::cerr << "^^^^^^^^^^" << std::endl;
BPVEC bfac;
mul.clear();
std::back_insert_iterator< BPVEC > bfac_bi( bfac );
mul_bi = std::back_insert_iterator<IVEC>(mul);
// std::cerr << CGAL::POLYNOMIAL::square_free_factorize( bp, bfac_bi, mul_bi ) << std::endl;
// std::cerr << bfac[0] << std::endl;
}
template <class FNT>
void integr() {
typedef CGAL::Polynomial<FNT> FPOLY;
typedef typename CGAL::Fraction_traits<FPOLY>::Numerator_type IPOLY;
typedef typename CGAL::Fraction_traits<FPOLY>::Denominator_type DENOM;
typename CGAL::Fraction_traits<FPOLY>::Decompose decompose;
typename CGAL::Fraction_traits<FPOLY>::Compose compose; (void) compose;
typedef typename IPOLY::NT INT;
// making polynomials integral and fractional
FPOLY fp(FNT(1)/FNT(2), FNT(2)/FNT(3), FNT(3)/FNT(5));
IPOLY ip; DENOM d;
decompose(fp, ip,d);
assert( d == DENOM(30) );
assert( ip == IPOLY(INT(15), INT(20), INT(18)) );
assert( fp == compose(ip, d) );
}
template <class NT>
void basic_tests() {
// tests that should work for all number types
typedef CGAL::Polynomial<NT> POLY;
datastru<NT>();
datastru<POLY>();
arithmetic<NT>();
arithmetic<POLY>();
io<NT>();
io<POLY>();
canon<NT>(typename CGAL::Algebraic_structure_traits<NT >::Algebraic_category());
#ifdef NiX_POLY_USE_NT_TESTS
// TODO: Replaced Is_real_comparable with RET::Is_real_embeddable, OK?
nt_tests<NT >(typename CGAL::Real_embeddable_traits<NT >::Is_real_embeddable());
// TODO: Replaced Is_real_comparable with RET::Is_real_embeddable, OK?
nt_tests<POLY>(typename CGAL::Real_embeddable_traits<POLY>::Is_real_embeddable());
#endif // NiX_POLY_USE_NT_TESTS
}
template <class NT>
void exact_tests() {
// tests that probably don't work for inexact or bounded number types
unigcdres<NT>(typename CGAL::Algebraic_structure_traits<NT>::Algebraic_category());
bigcdres<NT>(typename CGAL::Algebraic_structure_traits<NT>::Algebraic_category());
sqff<NT>();
}
template < class Rational >
void test_coefficients_to() {
typedef CGAL::Polynomial< int > Poly_i1;
typedef CGAL::Polynomial< Poly_i1 > Poly_i2;
typedef CGAL::Polynomial< Poly_i2 > Poly_i3;
typedef CGAL::Polynomial< Rational > Poly_rat1;
typedef CGAL::Polynomial< Poly_rat1 > Poly_rat2;
typedef CGAL::Polynomial< Poly_rat2 > Poly_rat3;
// univariate
{
Poly_i1 p = Poly_i1(-3,4,5,-9,0,2);
Poly_rat1 r = Poly_rat1(Rational(-3),Rational(4),Rational(5),
Rational(-9),Rational(0),Rational(2));
Poly_rat1 tmp;
CGAL::convert_to(p,tmp);
assert(tmp == r);
}
// bivariate
{
Poly_i1 p1 = Poly_i1( 3,-1, 2);
Poly_i1 p2 = Poly_i1(-5, 2,-4);
Poly_i2 p = Poly_i2(p1, p2);
Poly_rat1 r1 = Poly_rat1(Rational( 3),Rational(-1),Rational( 2));
Poly_rat1 r2 = Poly_rat1(Rational(-5),Rational( 2),Rational(-4));
Poly_rat2 r = Poly_rat2(r1,r2);
Poly_rat2 tmp;
CGAL::convert_to(p,tmp);
assert(tmp == r);
}
// trivariate
{
Poly_i1 p11 = Poly_i1( 7,-6, 9);
Poly_i1 p12 = Poly_i1(-1, 3);
Poly_i2 p1 = Poly_i2(p11, p12);
Poly_i1 p21 = Poly_i1(-17, 16, -19);
Poly_i1 p22 = Poly_i1( 11, -13);
Poly_i2 p2 = Poly_i2(p21, p22);
Poly_i3 p = Poly_i3(p1,p2);
Poly_rat1 r11 = Poly_rat1(Rational( 7),Rational(-6),Rational( 9));
Poly_rat1 r12 = Poly_rat1(Rational(-1),Rational( 3));
Poly_rat2 r1 = Poly_rat2(r11, r12);
Poly_rat1 r21 = Poly_rat1(Rational(-17),Rational( 16),Rational(-19));
Poly_rat1 r22 = Poly_rat1(Rational( 11),Rational(-13));
Poly_rat2 r2 = Poly_rat2(r21, r22);
Poly_rat3 r = Poly_rat3(r1,r2);
Poly_rat3 tmp;
CGAL::convert_to(p,tmp);
assert(tmp == r);
}
}
template <class AT>
void test_evaluate(){
CGAL_DEFINE_TYPES_FROM_AK(AT);
{
Poly_int1 P(3,0,2);
Integer x(2);
assert(P.evaluate(x)==Integer(11));
}{
Poly_rat1 P(3,0,2);
Integer x(2);
assert(P.evaluate(x)==Rational(11));
}{
Poly_int1 P(3,0,2);
Rational x(2);
assert(P.evaluate(x)==Rational(11));
}{
Poly_int1 P(3,0,2);
CGAL::Interval_nt<true> x(2);
assert(P.evaluate(x).inf()<=11);
assert(P.evaluate(x).sup()>=11);
}{
CGAL::Polynomial< CGAL::Interval_nt<true> > P(3,0,2);
CGAL::Interval_nt<true> x(2);
assert(P.evaluate(x).inf()<=11);
assert(P.evaluate(x).sup()>=11);
}{
CGAL::Polynomial< CGAL::Interval_nt<true> > P(3,0,2);
Integer x(2);
assert(P.evaluate(x).inf()<=11);
assert(P.evaluate(x).sup()>=11);
}
}
template <class AT>
void test_evaluate_homogeneous(){
CGAL_DEFINE_TYPES_FROM_AK(AT);
{
Poly_int1 P(3,0,2);
Integer u(2);
Integer v(1);
assert(P.evaluate_homogeneous(u,v)==Integer(11));
}{
Poly_rat1 P(2,-1,3);
Rational u(3);
Rational v(5);
assert(P.evaluate_homogeneous(u,v)==Rational(62));
assert(P.evaluate(u/v)==Rational(62)/Rational(25));
}
}
template <class NT>
void flat_iterator_tests() {
typedef CGAL::Polynomial<NT> POLY1;
typedef CGAL::Polynomial<POLY1> POLY2;
typedef CGAL::Polynomial<POLY2> POLY3;
POLY1
p1(NT( 1), NT( 2), NT( 3)),
p2(NT( 4), NT( 5), NT( 6)),
p3(NT( 7), NT( 8), NT( 9)),
p4(NT(10), NT(11), NT(12)),
p5(NT(13), NT(14), NT(15)),
p6(NT(16), NT(17), NT(18)),
p7(NT(19), NT(20), NT(21)),
p8(NT(22), NT(23), NT(24)),
p9(NT(25), NT(26), NT(27));
POLY2
q1(p1, p2, p3),
q2(p4, p5, p6),
q3(p7, p8, p9);
POLY3 r(q1, q2, q3);
typedef CGAL::Polynomial_traits_d<POLY1> PT1;
typedef CGAL::Polynomial_traits_d<POLY2> PT2;
typedef CGAL::Polynomial_traits_d<POLY3> PT3;
int i;
typename PT1::Innermost_coefficient_const_iterator it1; (void) it1;
typename PT1::Innermost_coefficient_const_begin begin1; (void) begin1;
typename PT1::Innermost_coefficient_const_end end1; (void) end1;
typename PT2::Innermost_coefficient_const_iterator it2; (void) it2;
typename PT2::Innermost_coefficient_const_begin begin2; (void) begin2;
typename PT2::Innermost_coefficient_const_end end2; (void) end2;
typename PT3::Innermost_coefficient_const_iterator it3; (void) it3;
typename PT3::Innermost_coefficient_const_begin begin3; (void) begin3;
typename PT3::Innermost_coefficient_const_end end3; (void) end3;
for (i = 1, it1 = begin1(p1); i <= 3; ++i, ++it1)
assert(*it1 == i);
assert(it1 == end1(p1));
for (i = 1, it2 = begin2(q1); i <= 9; ++i, ++it2)
assert(*it2 == i);
assert(it2 == end2(q1));
for (i = 1, it3 = begin3(r); i <= 27; ++i, ++it3)
assert(*it3 == i);
assert(it3 == end3(r));
}
void test_total_degree(){
typedef CGAL::Polynomial<int> Poly_1;
typedef CGAL::Polynomial<Poly_1> Poly_2;
assert(CGAL::total_degree(5) == 0);
assert(CGAL::total_degree(Poly_1(0)) == Poly_1(0).degree());
assert(CGAL::total_degree(Poly_1(1)) == Poly_1(1).degree());
assert(CGAL::total_degree(Poly_1(1,1)) == Poly_1(1,1).degree());
assert(CGAL::total_degree(Poly_1(0,0,0,1))== Poly_1(0,0,0,1).degree());
assert(CGAL::total_degree(Poly_2(0)) == 0);
assert(CGAL::total_degree(Poly_2(1)) == 0);
assert(CGAL::total_degree(Poly_2(Poly_1(1),Poly_1(1))) == 1);
assert(CGAL::total_degree(Poly_2(0,0,0,1))== 3);
assert(CGAL::total_degree(Poly_2(Poly_1(1,1),Poly_1(1))) == 1);
assert(CGAL::total_degree(Poly_2(0,0,1))== 2);
assert(CGAL::total_degree(Poly_2(Poly_1(0),Poly_1(0),Poly_1(1)))== 2);
assert(CGAL::total_degree(Poly_2(Poly_1(0),Poly_1(0),Poly_1(1,1)))== 3);
assert(CGAL::total_degree(Poly_2(Poly_1(1),Poly_1(0),Poly_1(1,1)))== 3);
assert(CGAL::total_degree(Poly_2(Poly_1(1),Poly_1(1),Poly_1(1,1)))== 3);
assert(CGAL::total_degree(Poly_2(Poly_1(1),Poly_1(1,1,1,1),Poly_1(1,1)))== 4);
}
template<class AT>
void test_scalar_factor_traits(){
{
typedef typename AT::Integer Integer;
typedef CGAL::Polynomial<Integer> Polynomial;
typedef CGAL::Scalar_factor_traits<Polynomial> SFT;
typedef typename AT::Integer Scalar;
typedef typename SFT::Scalar Scalar_;
BOOST_STATIC_ASSERT((::boost::is_same<Scalar_, Scalar>::value));
typename SFT::Scalar_factor sfac;
assert(sfac(Polynomial(0))==Scalar(0));
assert(sfac(Polynomial(9))==Scalar(9));
assert(sfac(Polynomial(9,15,30))==Scalar(3));
assert(sfac(Polynomial(0), Integer(0))==Scalar(0));
assert(sfac(Polynomial(0), Integer(1))==Scalar(1));
assert(sfac(Polynomial(0), Integer(2))==Scalar(2));
assert(sfac(Polynomial(9), Integer(0))==Scalar(9));
assert(sfac(Polynomial(9), Integer(6))==Scalar(3));
assert(sfac(Polynomial(15,0 ,30) , Integer(9))==Scalar(3));
assert(sfac(Polynomial(0 ,15,30) , Integer(9))==Scalar(3));
assert(sfac(Polynomial(18,15,30) , Integer(0))==Scalar(3));
}{
typedef typename AT::Integer Integer;
typedef CGAL::Sqrt_extension<Integer,Integer> EXT_1;
typedef CGAL::Polynomial<EXT_1 > Poly_1_ext_1;
typedef CGAL::Polynomial<Poly_1_ext_1> Poly_2_ext_1;
typedef CGAL::Scalar_factor_traits<Poly_2_ext_1> SFT;
typedef typename AT::Integer Scalar;
typedef typename SFT::Scalar Scalar_;
BOOST_STATIC_ASSERT((::boost::is_same<Scalar_, Scalar>::value));
typename SFT::Scalar_factor sfac;
assert(sfac(Poly_2_ext_1( ))==Integer(0));
assert(sfac(Poly_2_ext_1(1))==Integer(1));
assert(sfac(Poly_2_ext_1(2))==Integer(2));
EXT_1 a(18,27,456);
EXT_1 b( 0,15,456);
Poly_2_ext_1 p(Poly_1_ext_1(a,b),Poly_1_ext_1(b,a));
assert(sfac(p)==Integer(3));
}
}
template<class Polynomial_d, class T>
void test_interoperable_with(){
// construction
assert(Polynomial_d(-2) == T(-2));
// operators
Polynomial_d f(4);
assert( f + T(2) == Polynomial_d( 6));
assert( f - T(2) == Polynomial_d( 2));
assert( f * T(2) == Polynomial_d( 8));
assert( f / T(2) == Polynomial_d( 2));
assert( T(2) + f == Polynomial_d( 6));
assert( T(2) - f == Polynomial_d(-2));
assert( T(2) * f == Polynomial_d( 8));
assert( T(8) / f == Polynomial_d( 2));
f+=T(2);
f-=T(4);
f*=T(3);
f/=T(2);
assert( f == Polynomial_d(3));
assert( T(2) == Polynomial_d(2));
assert( T(2) != Polynomial_d(4));
assert( (T(2) < Polynomial_d(4)));
assert(!(T(4) < Polynomial_d(4)));
assert(!(T(6) < Polynomial_d(4)));
assert(!(T(2) > Polynomial_d(4)));
assert(!(T(4) > Polynomial_d(4)));
assert( (T(6) > Polynomial_d(4)));
assert( T(2) <= Polynomial_d(4));
assert( T(4) <= Polynomial_d(4));
assert(!(T(6) <= Polynomial_d(4)));
assert(!(T(2) >= Polynomial_d(4)));
assert( (T(4) >= Polynomial_d(4)));
assert( (T(6) >= Polynomial_d(4)));
}
template<class Polynomial_d>
void test_interoperable_poly(){
typedef CGAL::Polynomial_traits_d<Polynomial_d> PT;
typedef typename PT::Coefficient_type Coefficient_type;
typedef typename PT::Innermost_coefficient_type Innermost_coefficient_type;
test_interoperable_with<Polynomial_d,int>();
test_interoperable_with<Polynomial_d,Innermost_coefficient_type>();
test_interoperable_with<Polynomial_d,Coefficient_type>();
test_interoperable_with<Polynomial_d,Polynomial_d>();
}
template<class AT>
void test_interoperable_at(){
typedef typename AT::Integer Integer;
typedef CGAL::Sqrt_extension<Integer,Integer> EXT;
{
typedef int Coefficient_type;
typedef CGAL::Polynomial<Coefficient_type> Poly_1;
typedef CGAL::Polynomial<Poly_1> Poly_2;
typedef CGAL::Polynomial<Poly_2> Poly_3;
typedef CGAL::Polynomial<Poly_3> Poly_4;
test_interoperable_poly<Poly_1>();
test_interoperable_poly<Poly_2>();
// test_interoperable_poly<Poly_3>();
// test_interoperable_poly<Poly_4>();
}{
typedef Integer Coefficient_type;
typedef CGAL::Polynomial<Coefficient_type> Poly_1;
typedef CGAL::Polynomial<Poly_1> Poly_2;
typedef CGAL::Polynomial<Poly_2> Poly_3;
typedef CGAL::Polynomial<Poly_3> Poly_4;
test_interoperable_poly<Poly_1>();
test_interoperable_poly<Poly_2>();
// test_interoperable_poly<Poly_3>();
// test_interoperable_poly<Poly_4>();
}{
typedef EXT Coefficient_type;
typedef CGAL::Polynomial<Coefficient_type> Poly_1;
typedef CGAL::Polynomial<Poly_1> Poly_2;
typedef CGAL::Polynomial<Poly_2> Poly_3;
typedef CGAL::Polynomial<Poly_3> Poly_4;
test_interoperable_poly<Poly_1>();
test_interoperable_poly<Poly_2>();
// test_interoperable_poly<Poly_3>();
// test_interoperable_poly<Poly_4>();
}
}
template <class AT>
void test_AT(){
test_interoperable_at<AT>();
basic_tests<int>();
basic_tests<double>();
// TODO: The flat_iterator_test leads to compile errors. Need to be fixed
flat_iterator_tests<int>();
{
typedef typename AT::Integer Integer;
typedef typename CGAL::Polynomial<Integer> Polynomial;
typedef CGAL::Unique_factorization_domain_tag Tag;
typedef CGAL::Tag_true Is_exact;
CGAL::test_algebraic_structure<Polynomial,Tag, Is_exact>();
basic_tests<Integer>();
exact_tests<Integer>();
}{
typedef typename AT::Rational Rational;
typedef typename CGAL::Polynomial<Rational> Polynomial;
typedef CGAL::Euclidean_ring_tag Tag;
typedef CGAL::Tag_true Is_exact;
//can't use this test for Polynomials
//CGAL::test_algebraic_structure<Polynomial,Tag, Is_exact>();
basic_tests<Rational>();
exact_tests<Rational>();
integr<Rational>();
}
// TODO: This also leads to trouble with flattening. Needs to be fixed
psqff<AT>();
test_evaluate<AT>();
test_evaluate_homogeneous<AT>();
test_total_degree();
//test_scalar_factor_traits<AT>();
}
int main() {
CGAL::set_pretty_mode(std::cout);
// The MODULAR_TRAITS specializations for Polynomial
// =========================================================================
#ifdef CGAL_USE_LEDA
/* {
typedef CGAL::Polynomial<leda::integer> TESTT;
typedef CGAL::Modular MOD;
MOD::set_current_prime(7);
typedef CGAL::Modular_traits<TESTT> MOD_traits;
typedef MOD_traits::Residue_type MOD_NT;
MOD_traits::Modular_image MOD_image;
//definition of MOD::NT
assert((::boost::is_same<TESTT,
MOD_traits::NT>::value));
// Is_modularizable
assert(
(::boost::is_same<::CGAL::Tag_true,
MOD_traits::Is_modularizable>::value)
);
assert(MOD_image(TESTT(21)) == MOD_NT(0));
assert(MOD_image(TESTT(22)) == MOD_NT(1));
assert(MOD_image(TESTT(777777722)) == MOD_NT(1));
TESTT p(12,2,4,21),q(5,2,4);
assert(MOD_image(p)==MOD_image(q));
p=TESTT(7,21,14,-28);
assert(MOD_image(p)==MOD_NT(0));
}*/
// NT_Traits<POLY>::to_Interval
/* {
typedef leda_integer NT;
typedef CGAL::Interval Interval;
typedef CGAL::Polynomial<NT> POLY;
typedef CGAL::Polynomial<Interval> IPOLY;
typedef CGAL::Real_embeddable_traits<POLY>::To_interval::result_type RESULTTYPE;
CGAL::Real_embeddable_traits<POLY>::To_interval to_Interval;
BOOST_STATIC_ASSERT((::boost::is_same<IPOLY,RESULTTYPE>::value));
POLY p(NT(6), NT(-5), NT(1));
IPOLY ip=to_Interval(p);
assert(CGAL::in(12.0, ip.evaluate(Interval(-1))));
assert(CGAL::in(6.0 , ip.evaluate(Interval(0))));
assert(CGAL::in(2.0 , ip.evaluate(Interval(1))));
assert(CGAL::in(0.0 , ip.evaluate(Interval(2))));
assert(CGAL::in(0.0 , ip.evaluate(Interval(3))));
}
{ //POLY.evaluate(Interval)
typedef leda::integer NT;
typedef CGAL::Polynomial<NT> POLY;
typedef CGAL::Interval Interval;
POLY p(NT(6), NT(-5), NT(1));
assert(CGAL::in(12.0, p.evaluate(Interval(-1))));
assert(CGAL::in(6.0 , p.evaluate(Interval(0))));
assert(CGAL::in(2.0 , p.evaluate(Interval(1))));
assert(CGAL::in(0.0 , p.evaluate(Interval(2))));
assert(CGAL::in(0.0 , p.evaluate(Interval(3))));
assert(CGAL::in(12.0, p.evaluate_absolute(Interval(1))));
assert(CGAL::in(12.0, p.evaluate_absolute(Interval(-1))));
assert(CGAL::in(6.0 , p.evaluate_absolute(Interval(0))));
}*/
// coefficients_to
// test_coefficients_to<leda::rational>();
#endif // CGAL_USE_LEDA
#ifdef CGAL_USE_LEDA
{
typedef CGAL::LEDA_arithmetic_kernel AT;
test_AT<AT>();
}
#endif
#ifdef CGAL_USE_CORE
{
typedef CGAL::CORE_arithmetic_kernel AT;
test_AT<AT>();
}
#endif
return 0;
}
Reference in New Issue View Git Blame Copy Permalink