My Project
Public Types | Public Member Functions | Private Member Functions | Private Attributes
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_complexoperator[] (const int i)
 
gmp_complexevPointCoord (const int i)
 
gmp_complexgetRoot (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 ad0..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

Enumerator
none 
cspecial 
cspecialmu 
det 
onepoly 

Definition at line 68 of file mpr_numeric.h.

Constructor & Destructor Documentation

◆ rootContainer() [1/2]

rootContainer::rootContainer ( )

Definition at line 264 of file mpr_numeric.cc.

265 {
266  rt=none;
267 
268  coeffs= NULL;
269  ievpoint= NULL;
270  theroots= NULL;
271 
272  found_roots= false;
273 }
rootType rt
Definition: mpr_numeric.h:137
gmp_complex ** theroots
Definition: mpr_numeric.h:139
number * ievpoint
Definition: mpr_numeric.h:136
The main handler for Singular numbers which are suitable for Singular polynomials.
#define NULL
Definition: omList.c:12

◆ ~rootContainer()

rootContainer::~rootContainer ( )

Definition at line 277 of file mpr_numeric.cc.

278 {
279  int i;
280  // free coeffs, ievpoint
281  if ( ievpoint != NULL )
282  {
283  for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
284  omFreeSize( (void *)ievpoint, (anz+2) * sizeof( number ) );
285  }
286 
287  for ( i=0; i <= tdg; i++ )
288  if (coeffs[i]!=NULL) nDelete( coeffs + i );
289  omFreeSize( (void *)coeffs, (tdg+1) * sizeof( number ) );
290 
291  // theroots loeschen
292  for ( i=0; i < tdg; i++ ) delete theroots[i];
293  omFreeSize( (void *) theroots, (tdg)*sizeof(gmp_complex*) );
294 
295  //mprPROTnl("~rootContainer()");
296 }
int i
Definition: cfEzgcd.cc:132
gmp_complex numbers based on
Definition: mpr_complex.h:179
#define nDelete(n)
Definition: numbers.h:16
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260

◆ 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.

618 {
619  if(abs(x->imag())<abs(x->real())*e)
620  {
621  x->imag(0.0);
622  }
623 }
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
Variable x
Definition: cfModGcd.cc:4082

◆ 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.

804 {
805  int k;
806 
807  f0= *a[m];
808  ef= abs(f0);
809  f1= gmp_complex(0.0);
810  f2= f1;
811  ex= abs(x);
812 
813  for ( k= m-1; k >= 0; k-- )
814  {
815  f2 = ( x * f2 ) + f1;
816  f1 = ( x * f1 ) + f0;
817  f0 = ( x * f0 ) + *a[k];
818  ef = abs( f0 ) + ( ex * ef );
819  }
820 }
int m
Definition: cfEzgcd.cc:128
int k
Definition: cfEzgcd.cc:99

◆ 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.

825 {
826  int k;
827 
828  f0= *a[0];
829  ef= abs(f0);
830  f1= gmp_complex(0.0);
831  f2= f1;
832  ex= abs(x);
833 
834  for ( k= 1; k <= m; k++ )
835  {
836  f2 = ( x * f2 ) + f1;
837  f1 = ( x * f1 ) + f0;
838  f0 = ( x * f0 ) + *a[k];
839  ef = abs( f0 ) + ( ex * ef );
840  }
841 }

◆ divlin()

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

Definition at line 638 of file mpr_numeric.cc.

639 {
640  int i;
641  gmp_float o(1.0);
642 
643  if (abs(x)<o)
644  {
645  for (i= j-1; i > 0; i-- )
646  *a[i] += (*a[i+1]*x);
647  for (i= 0; i < j; i++ )
648  *a[i] = *a[i+1];
649  }
650  else
651  {
652  gmp_complex y(o/x);
653  for (i= 1; i < j; i++)
654  *a[i] += (*a[i-1]*y);
655  }
656 }
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
int j
Definition: facHensel.cc:110

◆ divquad()

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

Definition at line 658 of file mpr_numeric.cc.

659 {
660  int i;
661  gmp_float o(1.0),p(x.real()+x.real()),
662  q((x.real()*x.real())+(x.imag()*x.imag()));
663 
664  if (abs(x)<o)
665  {
666  *a[j-1] += (*a[j]*p);
667  for (i= j-2; i > 1; i-- )
668  *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
669  for (i= 0; i < j-1; i++ )
670  *a[i] = *a[i+2];
671  }
672  else
673  {
674  p = p/q;
675  q = o/q;
676  *a[1] += (*a[0]*p);
677  for (i= 2; i < j-1; i++)
678  *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
679  }
680 }
int p
Definition: cfModGcd.cc:4078

◆ evPointCoord()

gmp_complex & rootContainer::evPointCoord ( const int  i)

Definition at line 388 of file mpr_numeric.cc.

389 {
390  if (! ((i >= 0) && (i < anz+2) ) )
391  WarnS("rootContainer::evPointCoord: index out of range");
392  if (ievpoint == NULL)
393  WarnS("rootContainer::evPointCoord: ievpoint == NULL");
394 
395  if ( (rt == cspecialmu) && found_roots ) // FIX ME
396  {
397  if ( ievpoint[i] != NULL )
398  {
399  gmp_complex *tmp= new gmp_complex();
400  *tmp= numberToComplex(ievpoint[i], currRing->cf);
401  return *tmp;
402  }
403  else
404  {
405  Warn("rootContainer::evPointCoord: NULL index %d",i);
406  }
407  }
408 
409  // warning
410  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
411  gmp_complex *tmp= new gmp_complex();
412  return *tmp;
413 }
#define Warn
Definition: emacs.cc:77
#define WarnS
Definition: emacs.cc:78
gmp_complex numberToComplex(number num, const coeffs r)
Definition: mpr_complex.h:312
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13

◆ 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.

303 {
304  int i;
305  number nn= nInit(0);
306  var=_var;
307  tdg=_tdg;
308  coeffs=_coeffs;
309  rt=_rt;
310  anz=_anz;
311 
312  for ( i=0; i <= tdg; i++ )
313  {
314  if ( nEqual(coeffs[i],nn) )
315  {
316  nDelete( &coeffs[i] );
317  coeffs[i]=NULL;
318  }
319  }
320  nDelete( &nn );
321 
322  if ( rt == cspecialmu && _ievpoint ) // copy ievpoint
323  {
324  ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
325  for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
326  }
327 
328  theroots= NULL;
329  found_roots= false;
330 }
#define nEqual(n1, n2)
Definition: numbers.h:20
#define nCopy(n)
Definition: numbers.h:15
#define nInit(i)
Definition: numbers.h:24
#define omAlloc(size)
Definition: omAllocDecl.h:210

◆ getAnzElems()

int rootContainer::getAnzElems ( )
inline

Definition at line 95 of file mpr_numeric.h.

95 { return anz; }

◆ getAnzRoots()

int rootContainer::getAnzRoots ( )
inline

Definition at line 97 of file mpr_numeric.h.

97 { return tdg; }

◆ getLDim()

int rootContainer::getLDim ( )
inline

Definition at line 96 of file mpr_numeric.h.

96 { return anz; }

◆ getPoly()

poly rootContainer::getPoly ( )

Definition at line 334 of file mpr_numeric.cc.

335 {
336  int i;
337 
338  poly result= NULL;
339  poly ppos;
340 
341  if ( (rt == cspecial) || ( rt == cspecialmu ) )
342  {
343  for ( i= tdg; i >= 0; i-- )
344  {
345  if ( coeffs[i] )
346  {
347  poly p= pOne();
348  //pSetExp( p, var+1, i);
349  pSetExp( p, 1, i);
350  pSetCoeff( p, nCopy( coeffs[i] ) );
351  pSetm( p );
352  if (result)
353  {
354  ppos->next=p;
355  ppos=ppos->next;
356  }
357  else
358  {
359  result=p;
360  ppos=p;
361  }
362 
363  }
364  }
365  if (result!=NULL) pSetm( result );
366  }
367 
368  return result;
369 }
return result
Definition: facAbsBiFact.cc:75
#define pSetm(p)
Definition: polys.h:271
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pOne()
Definition: polys.h:315

◆ getRoot()

gmp_complex* rootContainer::getRoot ( const int  i)
inline

Definition at line 88 of file mpr_numeric.h.

89  {
90  return theroots[i];
91  }

◆ isfloat()

bool rootContainer::isfloat ( gmp_complex **  a)
private

Definition at line 625 of file mpr_numeric.cc.

626 {
627  gmp_float z(0.0);
628  gmp_complex *b;
629  for (int i=tdg; i >= 0; i-- )
630  {
631  b = &(*a[i]);
632  if (!(b->imag()==z))
633  return false;
634  }
635  return true;
636 }
CanonicalForm b
Definition: cfModGcd.cc:4103

◆ 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.

551 {
552  int iter,j;
553  gmp_float zero(0.0),one(1.0),deg(m);
554  gmp_float abx_g, err_g, fabs;
555  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
556  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
557 
558  gmp_float epss(0.1);
559  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),gmp_output_digits);
560 
561  for ( iter= 1; iter <= MAXIT; iter++ )
562  {
564  *its=iter;
565  if (type)
566  computefx(a,*x,m,b,d,f,abx_g,err_g);
567  else
568  computegx(a,*x,m,b,d,f,abx_g,err_g);
569  err_g *= epss; // EPSS;
570 
571  fabs = abs(b);
572  if (fabs <= err_g)
573  {
574  if ((fabs==zero) || (abs(d)==zero)) return;
575  *x -= (b/d); // a last newton-step
576  goto ende;
577  }
578 
579  g= d / b;
580  g2 = g * g;
581  h= g2 - (((f+f) / b ));
582  sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
583  gp= g + sq;
584  gm= g - sq;
585  if (abs(gp)<abs(gm))
586  {
587  dx = deg/gm;
588  }
589  else
590  {
591  if((gp.real()==zero)&&(gp.imag()==zero))
592  {
593  dx.real(cos((mprfloat)iter));
594  dx.imag(sin((mprfloat)iter));
595  dx = dx*(one+abx_g);
596  }
597  else
598  {
599  dx = deg/gp;
600  }
601  }
602  x1= *x - dx;
603 
604  if (*x == x1) goto ende;
605 
606  j = iter%MMOD;
607  if (j==0) j=MT;
608  if ( j % MT ) *x= x1;
609  else *x -= ( dx * frac_g[ j / MT ] );
610  }
611 
612  *its= MAXIT+1;
613 ende:
614  checkimag(x,epss);
615 }
g
Definition: cfModGcd.cc:4090
CanonicalForm gp
Definition: cfModGcd.cc:4102
FILE * f
Definition: checklibs.c:9
gmp_float imag() const
Definition: mpr_complex.h:235
gmp_float real() const
Definition: mpr_complex.h:234
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)
Definition: mpr_numeric.cc:822
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)
Definition: mpr_numeric.cc:801
void checkimag(gmp_complex *x, gmp_float &e)
Definition: mpr_numeric.cc:617
CFFListIterator iter
Definition: facAbsBiFact.cc:53
STATIC_VAR Poly * h
Definition: janet.cc:971
EXTERN_VAR size_t gmp_output_digits
Definition: mpr_base.h:115
gmp_float sin(const gmp_float &a)
Definition: mpr_complex.cc:333
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
gmp_float cos(const gmp_float &a)
Definition: mpr_complex.cc:338
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
#define ST_ROOTS_LG
Definition: mpr_global.h:82
double mprfloat
Definition: mpr_global.h:17
#define MT
Definition: mpr_numeric.cc:258
#define MMOD
Definition: mpr_numeric.cc:259
#define MR
Definition: mpr_numeric.cc:257
#define MAXIT
Definition: mpr_numeric.cc:260

◆ 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 ad0..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.

468 {
469  int i,j,k,its;
470  gmp_float zero(0.0);
471  gmp_complex x(0.0),o(1.0);
472  bool ret= true, isf=isfloat(a), type=true;
473 
474  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
475  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
476 
477  k = 0;
478  i = tdg;
479  j = i-1;
480  while (i>2)
481  {
482  // run laguer alg
483  x = zero;
484  laguer(ad, i, &x, &its, type);
485  if ( its > MAXIT )
486  {
487  type = !type;
488  x = zero;
489  laguer(ad, i, &x, &its, type);
490  }
491 
493  if ( its > MAXIT )
494  { // error
495  WarnS("Laguerre solver: Too many iterations!");
496  ret= false;
497  goto theend;
498  }
499  if ( polish )
500  {
501  laguer( a, tdg, &x, &its , type);
502  if ( its > MAXIT )
503  { // error
504  WarnS("Laguerre solver: Too many iterations in polish!");
505  ret= false;
506  goto theend;
507  }
508  }
509  if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x;
510  if (x.imag() == zero)
511  {
512  *roots[k] = x;
513  k++;
514  divlin(ad,x,i);
515  i--;
516  }
517  else
518  {
519  if(isf)
520  {
521  *roots[j] = x;
522  *roots[j-1]= gmp_complex(x.real(),-x.imag());
523  j -= 2;
524  divquad(ad,x,i);
525  i -= 2;
526  }
527  else
528  {
529  *roots[j] = x;
530  j--;
531  divlin(ad,x,i);
532  i--;
533  }
534  }
535  type = !type;
536  }
537  solvequad(ad,roots,k,j);
538  sortroots(roots,k,j,isf);
539 
540 theend:
541  mprSTICKYPROT("\n");
542  for ( i=0; i <= tdg; i++ ) delete ad[i];
543  omFreeSize( (void *) ad, (tdg+1)*sizeof( gmp_complex* ));
544 
545  return ret;
546 }
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 ...
Definition: mpr_numeric.cc:550
void divlin(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:638
void sortroots(gmp_complex **roots, int r, int c, bool isf)
Definition: mpr_numeric.cc:735
void divquad(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:658
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
Definition: mpr_numeric.cc:682
bool isfloat(gmp_complex **a)
Definition: mpr_numeric.cc:625
#define ST_ROOTS_LGSTEP
Definition: mpr_global.h:80

◆ operator[]()

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

Definition at line 82 of file mpr_numeric.h.

83  {
84  return *theroots[i];
85  }

◆ solvequad()

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

Definition at line 682 of file mpr_numeric.cc.

683 {
684  gmp_float zero(0.0);
685 
686  if ((j>k)
687  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
688  {
689  gmp_complex sq(zero);
690  gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
691  gmp_complex disk((h1 * h1) - h2);
692  if (disk.imag().isZero())
693  {
694  if (disk.real()<zero)
695  {
696  sq.real(zero);
697  sq.imag(sqrt(-disk.real()));
698  }
699  else
700  sq = (gmp_complex)sqrt(disk.real());
701  }
702  else
703  sq = sqrt(disk);
704  *r[k+1] = sq - h1;
705  sq += h1;
706  *r[k] = (gmp_complex)0.0-sq;
707  if(sq.imag().isZero())
708  {
709  k = j;
710  j++;
711  }
712  else
713  {
714  j = k;
715  k--;
716  }
717  }
718  else
719  {
720  if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero()))
721  {
722  WerrorS("precision lost, try again with higher precision");
723  }
724  else
725  {
726  *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
727  if(r[k]->imag().isZero())
728  j++;
729  else
730  k--;
731  }
732  }
733 }
bool isZero(const CFArray &A)
checks if entries of A are zero
void WerrorS(const char *s)
Definition: feFopen.cc:24

◆ solver()

bool rootContainer::solver ( const int  polishmode = PM_NONE)

Definition at line 437 of file mpr_numeric.cc.

438 {
439  int i;
440 
441  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
442  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
443  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
444 
445  // copy the coefficients of type number to type gmp_complex
446  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
447  for ( i=0; i <= tdg; i++ )
448  {
449  ad[i]= new gmp_complex();
450  if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i], currRing->cf );
451  }
452 
453  // now solve
454  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
455  if (!found_roots)
456  WarnS("rootContainer::solver: No roots found!");
457 
458  // free memory
459  for ( i=0; i <= tdg; i++ ) delete ad[i];
460  omFreeSize( (void *) ad, (tdg+1)*sizeof(gmp_complex*) );
461 
462  return found_roots;
463 }
bool laguer_driver(gmp_complex **a, gmp_complex **roots, bool polish=true)
Given the degree tdg and the tdg+1 complex coefficients ad0..tdg of the polynomial this routine succe...
Definition: mpr_numeric.cc:467

◆ sortre()

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

Definition at line 754 of file mpr_numeric.cc.

755 {
756  int pos,i;
757  gmp_complex *x,*y;
758 
759  pos = l;
760  x = r[pos];
761  for (i=l+inc; i<=u; i+=inc)
762  {
763  if (r[i]->real()<x->real())
764  {
765  pos = i;
766  x = r[pos];
767  }
768  }
769  if (pos>l)
770  {
771  if (inc==1)
772  {
773  for (i=pos; i>l; i--)
774  r[i] = r[i-1];
775  r[l] = x;
776  }
777  else
778  {
779  y = r[pos+1];
780  for (i=pos+1; i+1>l; i--)
781  r[i] = r[i-2];
782  if (x->imag()>y->imag())
783  {
784  r[l] = x;
785  r[l+1] = y;
786  }
787  else
788  {
789  r[l] = y;
790  r[l+1] = x;
791  }
792  }
793  }
794  else if ((inc==2)&&(x->imag()<r[l+1]->imag()))
795  {
796  r[l] = r[l+1];
797  r[l+1] = x;
798  }
799 }
int l
Definition: cfEzgcd.cc:100

◆ sortroots()

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

Definition at line 735 of file mpr_numeric.cc.

736 {
737  int j;
738 
739  for (j=0; j<r; j++) // the real roots
740  sortre(ro, j, r, 1);
741  if (c>=tdg) return;
742  if (isf)
743  {
744  for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly
745  sortre(ro, j, tdg-1, 2);
746  }
747  else
748  {
749  for (j=c; j+1<tdg; j++) // the complex roots for a general poly
750  sortre(ro, j, tdg-1, 1);
751  }
752 }
void sortre(gmp_complex **r, int l, int u, int inc)
Definition: mpr_numeric.cc:754

◆ swapRoots()

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

Definition at line 417 of file mpr_numeric.cc.

418 {
419  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
420  {
421  if ( to != from )
422  {
423  gmp_complex tmp( *theroots[from] );
424  *theroots[from]= *theroots[to];
425  *theroots[to]= tmp;
426  }
427  return true;
428  }
429 
430  // warning
431  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
432  return false;
433 }

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: