rootContainer Class Reference

complex root finder for univariate polynomials based on laguers algorithm More...

#include <mpr_numeric.h>

Public Types
enum	rootType {   none , cspecial , cspecialmu , det ,   onepoly }

Public Member Functions
	rootContainer ()
	~rootContainer ()
void	fillContainer (number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
bool	solver (const int polishmode=PM_NONE)
poly	getPoly ()
gmp_complex &	operator[] (const int i)
gmp_complex &	evPointCoord (const int i)
gmp_complex *	getRoot (const int i)
bool	swapRoots (const int from, const int to)
int	getAnzElems ()
int	getLDim ()
int	getAnzRoots ()

Private Member Functions
	rootContainer (const rootContainer &v)
bool	laguer_driver (gmp_complex **a, gmp_complex **roots, bool polish=true)
	Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg] (generated from the number coefficients coeffs[0..tdg]) of the polynomial this routine successively calls "laguer" and finds all m complex roots in roots[0..tdg]. More...
bool	isfloat (gmp_complex **a)
void	divlin (gmp_complex **a, gmp_complex x, int j)
void	divquad (gmp_complex **a, gmp_complex x, int j)
void	solvequad (gmp_complex **a, gmp_complex **r, int &k, int &j)
void	sortroots (gmp_complex **roots, int r, int c, bool isf)
void	sortre (gmp_complex **r, int l, int u, int inc)
void	laguer (gmp_complex **a, int m, gmp_complex *x, int *its, bool type)
	Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial, and given the complex value x, this routine improves x by Laguerre's method until it converges, within the achievable roundoff limit, to a root of the given polynomial. More...
void	computefx (gmp_complex **a, gmp_complex x, int m, gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, gmp_float &ex, gmp_float &ef)
void	computegx (gmp_complex **a, gmp_complex x, int m, gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, gmp_float &ex, gmp_float &ef)
void	checkimag (gmp_complex *x, gmp_float &e)

Private Attributes
int	var
int	tdg
number *	coeffs
number *	ievpoint
rootType	rt
gmp_complex **	theroots
int	anz
bool	found_roots

Detailed Description

complex root finder for univariate polynomials based on laguers algorithm

Definition at line 65 of file mpr_numeric.h.

Member Enumeration Documentation

rootType

enum rootContainer::rootType

Enumerator
none
cspecial
cspecialmu
det
onepoly

Definition at line 68 of file mpr_numeric.h.

{ none, cspecial, cspecialmu, det, onepoly };

Constructor & Destructor Documentation

rootContainer() [1/2]

rootContainer::rootContainer

(

)

Definition at line 264 of file mpr_numeric.cc.

{
  rt=none;
 
  coeffs= NULL;
  ievpoint= NULL;
  theroots= NULL;
 
  found_roots= false;
}

~rootContainer()

rootContainer::~rootContainer

(

)

Definition at line 277 of file mpr_numeric.cc.

{
  int i;
  // free coeffs, ievpoint
  if ( ievpoint != NULL )
  {
    for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
    omFreeSize( (void *)ievpoint, (anz+2) * sizeof( number ) );
  }
 
  for ( i=0; i <= tdg; i++ )
    if (coeffs[i]!=NULL) nDelete( coeffs + i );
  omFreeSize( (void *)coeffs, (tdg+1) * sizeof( number ) );
 
  // theroots loeschen
  for ( i=0; i < tdg; i++ ) delete theroots[i];
  omFreeSize( (void *) theroots, (tdg)*sizeof(gmp_complex*) );
 
  //mprPROTnl("~rootContainer()");
}

rootContainer() [2/2]

rootContainer::rootContainer ( const rootContainer &  v )

private

Member Function Documentation

checkimag()

void rootContainer::checkimag ( gmp_complex *  x, gmp_float &  e  )

private

Definition at line 617 of file mpr_numeric.cc.

{
  if(abs(x->imag())<abs(x->real())*e)
  {
    x->imag(0.0);
  }
}

computefx()

void rootContainer::computefx ( gmp_complex **  a, gmp_complex  x, int  m, gmp_complex &  f0, gmp_complex &  f1, gmp_complex &  f2, gmp_float &  ex, gmp_float &  ef  )

private

Definition at line 801 of file mpr_numeric.cc.

{
  int k;
 
  f0= *a[m];
  ef= abs(f0);
  f1= gmp_complex(0.0);
  f2= f1;
  ex= abs(x);
 
  for ( k= m-1; k >= 0; k-- )
  {
    f2 = ( x * f2 ) + f1;
    f1 = ( x * f1 ) + f0;
    f0 = ( x * f0 ) + *a[k];
    ef = abs( f0 ) + ( ex * ef );
  }
}

computegx()

void rootContainer::computegx ( gmp_complex **  a, gmp_complex  x, int  m, gmp_complex &  f0, gmp_complex &  f1, gmp_complex &  f2, gmp_float &  ex, gmp_float &  ef  )

private

Definition at line 822 of file mpr_numeric.cc.

{
  int k;
 
  f0= *a[0];
  ef= abs(f0);
  f1= gmp_complex(0.0);
  f2= f1;
  ex= abs(x);
 
  for ( k= 1; k <= m; k++ )
  {
    f2 = ( x * f2 ) + f1;
    f1 = ( x * f1 ) + f0;
    f0 = ( x * f0 ) + *a[k];
    ef = abs( f0 ) + ( ex * ef );
  }
}

divlin()

void rootContainer::divlin ( gmp_complex **  a, gmp_complex  x, int  j  )

private

Definition at line 638 of file mpr_numeric.cc.

{
  int i;
  gmp_float o(1.0);
 
  if (abs(x)<o)
  {
    for (i= j-1; i > 0; i-- )
      *a[i] += (*a[i+1]*x);
    for (i= 0; i < j; i++ )
      *a[i] = *a[i+1];
  }
  else
  {
    gmp_complex y(o/x);
    for (i= 1; i < j; i++)
      *a[i] += (*a[i-1]*y);
  }
}

divquad()

void rootContainer::divquad ( gmp_complex **  a, gmp_complex  x, int  j  )

private

Definition at line 658 of file mpr_numeric.cc.

{
  int i;
  gmp_float o(1.0),p(x.real()+x.real()),
            q((x.real()*x.real())+(x.imag()*x.imag()));
 
  if (abs(x)<o)
  {
    *a[j-1] += (*a[j]*p);
    for (i= j-2; i > 1; i-- )
      *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
    for (i= 0; i < j-1; i++ )
      *a[i] = *a[i+2];
  }
  else
  {
    p = p/q;
    q = o/q;
    *a[1] += (*a[0]*p);
    for (i= 2; i < j-1; i++)
      *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
  }
}

evPointCoord()

gmp_complex & rootContainer::evPointCoord

(

const int

i

)

Definition at line 388 of file mpr_numeric.cc.

{
  if (! ((i >= 0) && (i < anz+2) ) )
    WarnS("rootContainer::evPointCoord: index out of range");
  if (ievpoint == NULL)
    WarnS("rootContainer::evPointCoord: ievpoint == NULL");
 
  if ( (rt == cspecialmu) && found_roots ) // FIX ME
  {
    if ( ievpoint[i] != NULL )
    {
      gmp_complex *tmp= new gmp_complex();
      *tmp= numberToComplex(ievpoint[i], currRing->cf);
      return *tmp;
    }
    else
    {
      Warn("rootContainer::evPointCoord: NULL index %d",i);
    }
  }
 
  // warning
  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
  gmp_complex *tmp= new gmp_complex();
  return *tmp;
}

fillContainer()

void rootContainer::fillContainer	(	number *	_coeffs,
		number *	_ievpoint,
		const int	_var,
		const int	_tdg,
		const rootType	_rt,
		const int	_anz
	)

Definition at line 300 of file mpr_numeric.cc.

{
  int i;
  number nn= nInit(0);
  var=_var;
  tdg=_tdg;
  coeffs=_coeffs;
  rt=_rt;
  anz=_anz;
 
  for ( i=0; i <= tdg; i++ )
  {
    if ( nEqual(coeffs[i],nn) )
    {
      nDelete( &coeffs[i] );
      coeffs[i]=NULL;
    }
  }
  nDelete( &nn );
 
  if ( rt == cspecialmu && _ievpoint )  // copy ievpoint
  {
    ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
    for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
  }
 
  theroots= NULL;
  found_roots= false;
}

getAnzElems()

int rootContainer::getAnzElems ( )

inline

Definition at line 95 of file mpr_numeric.h.

{ return anz; }

getAnzRoots()

int rootContainer::getAnzRoots ( )

inline

Definition at line 97 of file mpr_numeric.h.

{ return tdg; }

getLDim()

int rootContainer::getLDim ( )

inline

Definition at line 96 of file mpr_numeric.h.

{ return anz; }

getPoly()

poly rootContainer::getPoly

(

)

Definition at line 334 of file mpr_numeric.cc.

{
  int i;
 
  poly result= NULL;
  poly ppos;
 
  if ( (rt == cspecial) || ( rt == cspecialmu ) )
  {
    for ( i= tdg; i >= 0; i-- )
    {
      if ( coeffs[i] )
      {
        poly p= pOne();
        //pSetExp( p, var+1, i);
        pSetExp( p, 1, i);
        pSetCoeff( p, nCopy( coeffs[i] ) );
        pSetm( p );
        if (result)
        {
          ppos->next=p;
          ppos=ppos->next;
        }
        else
        {
          result=p;
          ppos=p;
        }
 
      }
    }
    if (result!=NULL) pSetm( result );
  }
 
  return result;
}

getRoot()

gmp_complex * rootContainer::getRoot ( const int  i )

inline

Definition at line 88 of file mpr_numeric.h.

  {
    return theroots[i];
  }

isfloat()

bool rootContainer::isfloat ( gmp_complex **  a )

private

Definition at line 625 of file mpr_numeric.cc.

{
  gmp_float z(0.0);
  gmp_complex *b;
  for (int i=tdg; i >= 0; i-- )
  {
    b = &(*a[i]);
    if (!(b->imag()==z))
      return false;
  }
  return true;
}

laguer()

void rootContainer::laguer ( gmp_complex **  a, int  m, gmp_complex *  x, int *  its, bool  type  )

private

Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial, and given the complex value x, this routine improves x by Laguerre's method until it converges, within the achievable roundoff limit, to a root of the given polynomial.

The number of iterations taken is returned at its.

Definition at line 550 of file mpr_numeric.cc.

{
  int iter,j;
  gmp_float zero(0.0),one(1.0),deg(m);
  gmp_float abx_g, err_g, fabs;
  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
 
  gmp_float epss(0.1);
  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),gmp_output_digits);
 
  for ( iter= 1; iter <= MAXIT; iter++ )
  {
    mprSTICKYPROT(ST_ROOTS_LG);
    *its=iter;
    if (type)
      computefx(a,*x,m,b,d,f,abx_g,err_g);
    else
      computegx(a,*x,m,b,d,f,abx_g,err_g);
    err_g *= epss; // EPSS;
 
    fabs = abs(b);
    if (fabs <= err_g)
    {
      if ((fabs==zero) || (abs(d)==zero)) return;
      *x -= (b/d); // a last newton-step
      goto ende;
    }
 
    g= d / b;
    g2 = g * g;
    h= g2 - (((f+f) / b ));
    sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
    gp= g + sq;
    gm= g - sq;
    if (abs(gp)<abs(gm))
    {
      dx = deg/gm;
    }
    else
    {
      if((gp.real()==zero)&&(gp.imag()==zero))
      {
        dx.real(cos((mprfloat)iter));
        dx.imag(sin((mprfloat)iter));
        dx = dx*(one+abx_g);
      }
      else
      {
        dx = deg/gp;
      }
    }
    x1= *x - dx;
 
    if (*x == x1) goto ende;
 
    j = iter%MMOD;
    if (j==0) j=MT;
    if ( j % MT ) *x= x1;
    else *x -= ( dx * frac_g[ j / MT ] );
  }
 
  *its= MAXIT+1;
ende:
  checkimag(x,epss);
}

laguer_driver()

bool rootContainer::laguer_driver ( gmp_complex **  a, gmp_complex **  roots, bool  polish = true  )

private

Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg] (generated from the number coefficients coeffs[0..tdg]) of the polynomial this routine successively calls "laguer" and finds all m complex roots in roots[0..tdg].

The bool var "polish" should be input as "true" if polishing (also by "laguer") is desired, "false" if the roots will be subsequently polished by other means.

Definition at line 467 of file mpr_numeric.cc.

{
  int i,j,k,its;
  gmp_float zero(0.0);
  gmp_complex x(0.0),o(1.0);
  bool ret= true, isf=isfloat(a), type=true;
 
  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
 
  k = 0;
  i = tdg;
  j = i-1;
  while (i>2)
  {
    // run laguer alg
    x = zero;
    laguer(ad, i, &x, &its, type);
    if ( its > MAXIT )
    {
      type = !type;
      x = zero;
      laguer(ad, i, &x, &its, type);
    }
 
    mprSTICKYPROT(ST_ROOTS_LGSTEP);
    if ( its > MAXIT )
    {  // error
      WarnS("Laguerre solver: Too many iterations!");
      ret= false;
      goto theend;
    }
    if ( polish )
    {
      laguer( a, tdg, &x, &its , type);
      if ( its > MAXIT )
      {  // error
        WarnS("Laguerre solver: Too many iterations in polish!");
        ret= false;
        goto theend;
      }
    }
    if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x;
    if (x.imag() == zero)
    {
      *roots[k] = x;
      k++;
      divlin(ad,x,i);
      i--;
    }
    else
    {
      if(isf)
      {
        *roots[j] = x;
        *roots[j-1]= gmp_complex(x.real(),-x.imag());
        j -= 2;
        divquad(ad,x,i);
        i -= 2;
      }
      else
      {
        *roots[j] = x;
        j--;
        divlin(ad,x,i);
        i--;
      }
    }
    type = !type;
  }
  solvequad(ad,roots,k,j);
  sortroots(roots,k,j,isf);
 
theend:
  mprSTICKYPROT("\n");
  for ( i=0; i <= tdg; i++ ) delete ad[i];
  omFreeSize( (void *) ad, (tdg+1)*sizeof( gmp_complex* ));
 
  return ret;
}

operator[]()

gmp_complex & rootContainer::operator[] ( const int  i )

inline

Definition at line 82 of file mpr_numeric.h.

  {
    return *theroots[i];
  }

solvequad()

void rootContainer::solvequad ( gmp_complex **  a, gmp_complex **  r, int &  k, int &  j  )

private

Definition at line 682 of file mpr_numeric.cc.

{
  gmp_float zero(0.0);
 
  if ((j>k)
  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
  {
    gmp_complex sq(zero);
    gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
    gmp_complex disk((h1 * h1) - h2);
    if (disk.imag().isZero())
    {
      if (disk.real()<zero)
      {
        sq.real(zero);
        sq.imag(sqrt(-disk.real()));
      }
      else
        sq = (gmp_complex)sqrt(disk.real());
    }
    else
      sq = sqrt(disk);
    *r[k+1] = sq - h1;
    sq += h1;
    *r[k] = (gmp_complex)0.0-sq;
    if(sq.imag().isZero())
    {
      k = j;
      j++;
    }
    else
    {
      j = k;
      k--;
    }
  }
  else
  {
    if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero()))
    {
      WerrorS("precision lost, try again with higher precision");
    }
    else
    {
      *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
      if(r[k]->imag().isZero())
        j++;
      else
        k--;
    }
  }
}

solver()

bool rootContainer::solver

(

const int

polishmode = PM_NONE

)

Definition at line 437 of file mpr_numeric.cc.

{
  int i;
 
  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
 
  // copy the coefficients of type number to type gmp_complex
  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
  for ( i=0; i <= tdg; i++ )
  {
    ad[i]= new gmp_complex();
    if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i], currRing->cf );
  }
 
  // now solve
  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
  if (!found_roots)
    WarnS("rootContainer::solver: No roots found!");
 
  // free memory
  for ( i=0; i <= tdg; i++ ) delete ad[i];
  omFreeSize( (void *) ad, (tdg+1)*sizeof(gmp_complex*) );
 
  return found_roots;
}

sortre()

void rootContainer::sortre ( gmp_complex **  r, int  l, int  u, int  inc  )

private

Definition at line 754 of file mpr_numeric.cc.

{
  int pos,i;
  gmp_complex *x,*y;
 
  pos = l;
  x = r[pos];
  for (i=l+inc; i<=u; i+=inc)
  {
    if (r[i]->real()<x->real())
    {
      pos = i;
      x = r[pos];
    }
  }
  if (pos>l)
  {
    if (inc==1)
    {
      for (i=pos; i>l; i--)
        r[i] = r[i-1];
      r[l] = x;
    }
    else
    {
      y = r[pos+1];
      for (i=pos+1; i+1>l; i--)
        r[i] = r[i-2];
      if (x->imag()>y->imag())
      {
        r[l] = x;
        r[l+1] = y;
      }
      else
      {
        r[l] = y;
        r[l+1] = x;
      }
    }
  }
  else if ((inc==2)&&(x->imag()<r[l+1]->imag()))
  {
    r[l] = r[l+1];
    r[l+1] = x;
  }
}

sortroots()

void rootContainer::sortroots ( gmp_complex **  roots, int  r, int  c, bool  isf  )

private

Definition at line 735 of file mpr_numeric.cc.

{
  int j;
 
  for (j=0; j<r; j++) // the real roots
    sortre(ro, j, r, 1);
  if (c>=tdg) return;
  if (isf)
  {
    for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly
      sortre(ro, j, tdg-1, 2);
  }
  else
  {
    for (j=c; j+1<tdg; j++) // the complex roots for a general poly
      sortre(ro, j, tdg-1, 1);
  }
}

swapRoots()

bool rootContainer::swapRoots	(	const int	from,
		const int	to
	)

Definition at line 417 of file mpr_numeric.cc.

{
  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
  {
    if ( to != from )
    {
      gmp_complex tmp( *theroots[from] );
      *theroots[from]= *theroots[to];
      *theroots[to]= tmp;
    }
    return true;
  }
 
  // warning
  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
  return false;
}

Field Documentation

anz

int rootContainer::anz

private

Definition at line 141 of file mpr_numeric.h.

coeffs

number* rootContainer::coeffs

private

Definition at line 135 of file mpr_numeric.h.

found_roots

bool rootContainer::found_roots

private

Definition at line 142 of file mpr_numeric.h.

ievpoint

number* rootContainer::ievpoint

private

Definition at line 136 of file mpr_numeric.h.

rt

rootType rootContainer::rt

private

Definition at line 137 of file mpr_numeric.h.

tdg

int rootContainer::tdg

private

Definition at line 133 of file mpr_numeric.h.

theroots

gmp_complex** rootContainer::theroots

private

Definition at line 139 of file mpr_numeric.h.

var

int rootContainer::var

private

Definition at line 132 of file mpr_numeric.h.

---

The documentation for this class was generated from the following files:

• kernel/numeric/mpr_numeric.h
• kernel/numeric/mpr_numeric.cc