CGAL users discussion list

[Eric Berberich <>]

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