Skip to Content.
Sympa Menu

cgal-discuss - Re: [cgal-discuss] Possible convergence bug in CGAL algebraic kernel

Subject: CGAL users discussion list

List archive

Re: [cgal-discuss] Possible convergence bug in CGAL algebraic kernel


Chronological Thread 
  • From: Eric Berberich <>
  • To:
  • Cc: adel ahmadyan <>
  • Subject: Re: [cgal-discuss] Possible convergence bug in CGAL algebraic kernel
  • Date: Thu, 12 Jul 2012 08:15:16 +0200

Hi,

On 7/11/12 7:38 PM, adel ahmadyan wrote:

Hi,
I ran into a converge problem while using CGAL for solving polynomials.
Using algebraic kernel, CGAL says that the following poly does not have
any roots,
but it has 3 roots. It seems that CGAL does not converge on these roots
due to numerical problems. (CGAL's example for simpler polys are working
fine).

I can only guess of which type the coefficient are, as you have not included these lines.

Please, send files that we can execute from the hook, i.e. with ALL include and with all variable specified!

So I guess you you used:


#include <CGAL/config.h>
#include <CGAl/Algebraic_kernel_d_1.h>
#include <iostream>
using namespace std;

int main() {

double a_0 = -1.37937;
double a_1 = -0.384008;
double a_2 = 0.360753;
double a_3 = -0.020945;

// your posted code here

}


It indeed compiles and reports no roots.

The problem is not in the AK_1, it's in type conversion from double (a floating-point number type) to Gmpz (an integral number). If you use double coefficients, as I have to assume, your coefficients get *rounded* to the nearest integer.

Gmpz ca_0(a_0); // ca_0 is now -1
Gmpz ca_1(a_1); // ca_1 is now 0
Gmpz ca_2(a_2); // ca_2 is now 0
Gmpz ca_3(a_3); // ca_3 is now 0

I've added the following output:

cout << "poly: " << a_3*x*x*x + a_2*x*x+ a_1*x + a_0 << std::endl;

which shows:

poly: P[0(0,-1)]

i.e. you're trying next to find the real roots of a constant polynomial, which has none.


Beyond: Using double to represent exact (decimal) coefficient is a bad idea, because of the rounding errors that occur (I mean your code even contains a BIG rounding error). The correct way (not to loose any root) is storing the coefficient in an exact Rational number type and then constructing from them the integral polynomial needed (i.e. one with integer coefficient). Note that multiplying with a common denominator does not change the roots of a polynomial.

A quick hack for you is multiplying the coefficient manually with 10ˆ6.

int a_0 = -1379370;
int a_1 = -384008;
int a_2 = 360753;
int a_3 = -20945;

then it's correct

poly: P[3(0,-1379370)(1,-384008)(2,360753)(3,-20945)]
Number of roots is : 3
-1.45014 1
2.87438 1
15.7996 1



HTH,
eriC

poly = a0 + a1x + a2x^2 + a3x^3
a0 = -1.37937
a1 = -0.384008
a2 = 0.360753
a3 = -0.020945

This is CGAL/c++ code:

CGAL::Algebraic_kernel_d_1<CGAL::Gmpz> ak; // an object of
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Solve_1 solve_1 =
ak.solve_1_object();
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Polynomial_1 x =

CGAL::shift(CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Polynomial_1(1),1);
// the monomial x
std::vector<std::pair<
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Algebraic_real_1,
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Multiplicity_type>
> mroots;
cout << "a0= " << a_0 << endl << " a1= " << a_1 << endl <<
" a2= " << a_2 << endl << " a3=" << a_3 << endl ;
solve_1( a_3*x*x*x + a_2*x*x+ a_1*x + a_0,
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Bound(-100),
CGAL::Algebraic_kernel_d_1<CGAL::Gmpz>::Bound(100),
std::back_inserter(mroots));
std::cout << "Number of roots is : " <<
mroots.size() << "\n";
for(int i=0;i<mroots.size();i++){
cout << CGAL::to_double(mroots[i].first) << " " <<
CGAL::to_double(mroots[i].second) << endl ;
}

MATLAB:
p = [-0.020945 0.360753 -0.384008 -1.37937] ;
r = roots(p)
r =
15.7996
2.8744
-1.4501





Archive powered by MHonArc 2.6.18.

Top of Page