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

Enumerator
none 
cspecial 
cspecialmu 
det 
onepoly 

Definition at line 68 of file mpr_numeric.h.

Constructor & Destructor Documentation

rootContainer::rootContainer ( )

Definition at line 278 of file mpr_numeric.cc.

279 {
280  rt=none;
281 
282  coeffs= NULL;
283  ievpoint= NULL;
284  theroots= NULL;
285 
286  found_roots= false;
287 }
number * ievpoint
Definition: mpr_numeric.h:136
rootType rt
Definition: mpr_numeric.h:137
The main handler for Singular numbers which are suitable for Singular polynomials.
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define NULL
Definition: omList.c:10
rootContainer::~rootContainer ( )

Definition at line 291 of file mpr_numeric.cc.

292 {
293  int i;
294  // free coeffs, ievpoint
295  if ( ievpoint != NULL )
296  {
297  for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
298  omFreeSize( (void *)ievpoint, (anz+2) * sizeof( number ) );
299  }
300 
301  for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i );
302  omFreeSize( (void *)coeffs, (tdg+1) * sizeof( number ) );
303 
304  // theroots löschen
305  for ( i=0; i < tdg; i++ ) delete theroots[i];
306  omFreeSize( (void *) theroots, (tdg)*sizeof(gmp_complex*) );
307 
308  //mprPROTnl("~rootContainer()");
309 }
number * ievpoint
Definition: mpr_numeric.h:136
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
gmp_complex numbers based on
Definition: mpr_complex.h:178
The main handler for Singular numbers which are suitable for Singular polynomials.
int i
Definition: cfEzgcd.cc:123
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define nDelete(n)
Definition: numbers.h:16
#define NULL
Definition: omList.c:10
rootContainer::rootContainer ( const rootContainer v)
private

Member Function Documentation

void rootContainer::checkimag ( gmp_complex x,
gmp_float e 
)
private

Definition at line 630 of file mpr_numeric.cc.

631 {
632  if(abs(x->imag())<abs(x->real())*e)
633  {
634  x->imag(0.0);
635  }
636 }
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
gmp_float imag() const
Definition: mpr_complex.h:235
gmp_float real() const
Definition: mpr_complex.h:234
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 814 of file mpr_numeric.cc.

817 {
818  int k;
819 
820  f0= *a[m];
821  ef= abs(f0);
822  f1= gmp_complex(0.0);
823  f2= f1;
824  ex= abs(x);
825 
826  for ( k= m-1; k >= 0; k-- )
827  {
828  f2 = ( x * f2 ) + f1;
829  f1 = ( x * f1 ) + f0;
830  f0 = ( x * f0 ) + *a[k];
831  ef = abs( f0 ) + ( ex * ef );
832  }
833 }
int k
Definition: cfEzgcd.cc:93
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
int m
Definition: cfEzgcd.cc:119
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 835 of file mpr_numeric.cc.

838 {
839  int k;
840 
841  f0= *a[0];
842  ef= abs(f0);
843  f1= gmp_complex(0.0);
844  f2= f1;
845  ex= abs(x);
846 
847  for ( k= 1; k <= m; k++ )
848  {
849  f2 = ( x * f2 ) + f1;
850  f1 = ( x * f1 ) + f0;
851  f0 = ( x * f0 ) + *a[k];
852  ef = abs( f0 ) + ( ex * ef );
853  }
854 }
int k
Definition: cfEzgcd.cc:93
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
int m
Definition: cfEzgcd.cc:119
void rootContainer::divlin ( gmp_complex **  a,
gmp_complex  x,
int  j 
)
private

Definition at line 651 of file mpr_numeric.cc.

652 {
653  int i;
654  gmp_float o(1.0);
655 
656  if (abs(x)<o)
657  {
658  for (i= j-1; i > 0; i-- )
659  *a[i] += (*a[i+1]*x);
660  for (i= 0; i < j; i++ )
661  *a[i] = *a[i+1];
662  }
663  else
664  {
665  gmp_complex y(o/x);
666  for (i= 1; i < j; i++)
667  *a[i] += (*a[i-1]*y);
668  }
669 }
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
void rootContainer::divquad ( gmp_complex **  a,
gmp_complex  x,
int  j 
)
private

Definition at line 671 of file mpr_numeric.cc.

672 {
673  int i;
674  gmp_float o(1.0),p(x.real()+x.real()),
675  q((x.real()*x.real())+(x.imag()*x.imag()));
676 
677  if (abs(x)<o)
678  {
679  *a[j-1] += (*a[j]*p);
680  for (i= j-2; i > 1; i-- )
681  *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
682  for (i= 0; i < j-1; i++ )
683  *a[i] = *a[i+2];
684  }
685  else
686  {
687  p = p/q;
688  q = o/q;
689  *a[1] += (*a[0]*p);
690  for (i= 2; i < j-1; i++)
691  *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
692  }
693 }
return P p
Definition: myNF.cc:203
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
gmp_float imag() const
Definition: mpr_complex.h:235
gmp_float real() const
Definition: mpr_complex.h:234
gmp_complex & rootContainer::evPointCoord ( const int  i)

Definition at line 401 of file mpr_numeric.cc.

402 {
403  if (! ((i >= 0) && (i < anz+2) ) )
404  WarnS("rootContainer::evPointCoord: index out of range");
405  if (ievpoint == NULL)
406  WarnS("rootContainer::evPointCoord: ievpoint == NULL");
407 
408  if ( (rt == cspecialmu) && found_roots ) // FIX ME
409  {
410  if ( ievpoint[i] != NULL )
411  {
412  gmp_complex *tmp= new gmp_complex();
413  *tmp= numberToComplex(ievpoint[i], currRing->cf);
414  return *tmp;
415  }
416  else
417  {
418  Warn("rootContainer::evPointCoord: NULL index %d",i);
419  }
420  }
421 
422  // warning
423  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
424  gmp_complex *tmp= new gmp_complex();
425  return *tmp;
426 }
number * ievpoint
Definition: mpr_numeric.h:136
rootType rt
Definition: mpr_numeric.h:137
gmp_complex numbers based on
Definition: mpr_complex.h:178
#define WarnS
Definition: emacs.cc:81
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
gmp_complex numberToComplex(number num, const coeffs r)
Definition: mpr_complex.h:312
int i
Definition: cfEzgcd.cc:123
#define NULL
Definition: omList.c:10
#define Warn
Definition: emacs.cc:80
void rootContainer::fillContainer ( number *  _coeffs,
number *  _ievpoint,
const int  _var,
const int  _tdg,
const rootType  _rt,
const int  _anz 
)

Definition at line 313 of file mpr_numeric.cc.

316 {
317  int i;
318  number nn= nInit(0);
319  var=_var;
320  tdg=_tdg;
321  coeffs=_coeffs;
322  rt=_rt;
323  anz=_anz;
324 
325  for ( i=0; i <= tdg; i++ )
326  {
327  if ( nEqual(coeffs[i],nn) )
328  {
329  nDelete( &coeffs[i] );
330  coeffs[i]=NULL;
331  }
332  }
333  nDelete( &nn );
334 
335  if ( rt == cspecialmu && _ievpoint ) // copy ievpoint
336  {
337  ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
338  for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
339  }
340 
341  theroots= NULL;
342  found_roots= false;
343 }
number * ievpoint
Definition: mpr_numeric.h:136
rootType rt
Definition: mpr_numeric.h:137
#define nEqual(n1, n2)
Definition: numbers.h:20
#define omAlloc(size)
Definition: omAllocDecl.h:210
The main handler for Singular numbers which are suitable for Singular polynomials.
int i
Definition: cfEzgcd.cc:123
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define nDelete(n)
Definition: numbers.h:16
#define NULL
Definition: omList.c:10
#define nCopy(n)
Definition: numbers.h:15
#define nInit(i)
Definition: numbers.h:24
int rootContainer::getAnzElems ( )
inline

Definition at line 95 of file mpr_numeric.h.

95 { return anz; }
int rootContainer::getAnzRoots ( )
inline

Definition at line 97 of file mpr_numeric.h.

97 { return tdg; }
int rootContainer::getLDim ( )
inline

Definition at line 96 of file mpr_numeric.h.

96 { return anz; }
poly rootContainer::getPoly ( )

Definition at line 347 of file mpr_numeric.cc.

348 {
349  int i;
350 
351  poly result= NULL;
352  poly ppos;
353 
354  if ( (rt == cspecial) || ( rt == cspecialmu ) )
355  {
356  for ( i= tdg; i >= 0; i-- )
357  {
358  if ( coeffs[i] )
359  {
360  poly p= pOne();
361  //pSetExp( p, var+1, i);
362  pSetExp( p, 1, i);
363  pSetCoeff( p, nCopy( coeffs[i] ) );
364  pSetm( p );
365  if (result)
366  {
367  ppos->next=p;
368  ppos=ppos->next;
369  }
370  else
371  {
372  result=p;
373  ppos=p;
374  }
375 
376  }
377  }
378  if (result!=NULL) pSetm( result );
379  }
380 
381  return result;
382 }
#define pSetm(p)
Definition: polys.h:241
#define pSetExp(p, i, v)
Definition: polys.h:42
return P p
Definition: myNF.cc:203
rootType rt
Definition: mpr_numeric.h:137
The main handler for Singular numbers which are suitable for Singular polynomials.
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:286
#define NULL
Definition: omList.c:10
#define nCopy(n)
Definition: numbers.h:15
polyrec * poly
Definition: hilb.h:10
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
return result
Definition: facAbsBiFact.cc:76
gmp_complex* rootContainer::getRoot ( const int  i)
inline

Definition at line 88 of file mpr_numeric.h.

89  {
90  return theroots[i];
91  }
int i
Definition: cfEzgcd.cc:123
gmp_complex ** theroots
Definition: mpr_numeric.h:139
bool rootContainer::isfloat ( gmp_complex **  a)
private

Definition at line 638 of file mpr_numeric.cc.

639 {
640  gmp_float z(0.0);
641  gmp_complex *b;
642  for (int i=tdg; i >= 0; i-- )
643  {
644  b = &(*a[i]);
645  if (!(b->imag()==z))
646  return false;
647  }
648  return true;
649 }
gmp_complex numbers based on
Definition: mpr_complex.h:178
int i
Definition: cfEzgcd.cc:123
gmp_float imag() const
Definition: mpr_complex.h:235
const poly b
Definition: syzextra.cc:213
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 563 of file mpr_numeric.cc.

564 {
565  int iter,j;
566  gmp_float zero(0.0),one(1.0),deg(m);
567  gmp_float abx_g, err_g, fabs;
568  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
569  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
570 
571  gmp_float epss(0.1);
572  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),gmp_output_digits);
573 
574  for ( iter= 1; iter <= MAXIT; iter++ )
575  {
577  *its=iter;
578  if (type)
579  computefx(a,*x,m,b,d,f,abx_g,err_g);
580  else
581  computegx(a,*x,m,b,d,f,abx_g,err_g);
582  err_g *= epss; // EPSS;
583 
584  fabs = abs(b);
585  if (fabs <= err_g)
586  {
587  if ((fabs==zero) || (abs(d)==zero)) return;
588  *x -= (b/d); // a last newton-step
589  goto ende;
590  }
591 
592  g= d / b;
593  g2 = g * g;
594  h= g2 - (((f+f) / b ));
595  sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
596  gp= g + sq;
597  gm= g - sq;
598  if (abs(gp)<abs(gm))
599  {
600  dx = deg/gm;
601  }
602  else
603  {
604  if((gp.real()==zero)&&(gp.imag()==zero))
605  {
606  dx.real(cos((mprfloat)iter));
607  dx.imag(sin((mprfloat)iter));
608  dx = dx*(one+abx_g);
609  }
610  else
611  {
612  dx = deg/gp;
613  }
614  }
615  x1= *x - dx;
616 
617  if (*x == x1) goto ende;
618 
619  j = iter%MMOD;
620  if (j==0) j=MT;
621  if ( j % MT ) *x= x1;
622  else *x -= ( dx * frac_g[ j / MT ] );
623  }
624 
625  *its= MAXIT+1;
626 ende:
627  checkimag(x,epss);
628 }
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
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:835
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:814
#define MT
Definition: mpr_numeric.cc:272
f
Definition: cfModGcd.cc:4022
#define MAXIT
Definition: mpr_numeric.cc:274
CFFListIterator iter
Definition: facAbsBiFact.cc:54
void checkimag(gmp_complex *x, gmp_float &e)
Definition: mpr_numeric.cc:630
#define MMOD
Definition: mpr_numeric.cc:273
g
Definition: cfModGcd.cc:4031
gmp_complex numbers based on
Definition: mpr_complex.h:178
double mprfloat
Definition: mpr_global.h:17
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
int j
Definition: myNF.cc:70
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:329
int m
Definition: cfEzgcd.cc:119
#define ST_ROOTS_LG
Definition: mpr_global.h:82
gmp_float cos(const gmp_float &a)
Definition: mpr_complex.cc:340
CanonicalForm gp
Definition: cfModGcd.cc:4043
gmp_float imag() const
Definition: mpr_complex.h:235
#define MR
Definition: mpr_numeric.cc:271
size_t gmp_output_digits
Definition: mpr_complex.cc:44
static Poly * h
Definition: janet.cc:978
const poly b
Definition: syzextra.cc:213
gmp_float sin(const gmp_float &a)
Definition: mpr_complex.cc:335
gmp_float real() const
Definition: mpr_complex.h:234
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 480 of file mpr_numeric.cc.

481 {
482  int i,j,k,its;
483  gmp_float zero(0.0);
484  gmp_complex x(0.0),o(1.0);
485  bool ret= true, isf=isfloat(a), type=true;
486 
487  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
488  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
489 
490  k = 0;
491  i = tdg;
492  j = i-1;
493  while (i>2)
494  {
495  // run laguer alg
496  x = zero;
497  laguer(ad, i, &x, &its, type);
498  if ( its > MAXIT )
499  {
500  type = !type;
501  x = zero;
502  laguer(ad, i, &x, &its, type);
503  }
504 
506  if ( its > MAXIT )
507  { // error
508  WarnS("Laguerre solver: Too many iterations!");
509  ret= false;
510  goto theend;
511  }
512  if ( polish )
513  {
514  laguer( a, tdg, &x, &its , type);
515  if ( its > MAXIT )
516  { // error
517  WarnS("Laguerre solver: Too many iterations in polish!");
518  ret= false;
519  goto theend;
520  }
521  }
522  if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x;
523  if (x.imag() == zero)
524  {
525  *roots[k] = x;
526  k++;
527  divlin(ad,x,i);
528  i--;
529  }
530  else
531  {
532  if(isf)
533  {
534  *roots[j] = x;
535  *roots[j-1]= gmp_complex(x.real(),-x.imag());
536  j -= 2;
537  divquad(ad,x,i);
538  i -= 2;
539  }
540  else
541  {
542  *roots[j] = x;
543  j--;
544  divlin(ad,x,i);
545  i--;
546  }
547  }
548  type = !type;
549  }
550  solvequad(ad,roots,k,j);
551  sortroots(roots,k,j,isf);
552 
553 theend:
554  mprSTICKYPROT("\n");
555  for ( i=0; i <= tdg; i++ ) delete ad[i];
556  omFreeSize( (void *) ad, (tdg+1)*sizeof( gmp_complex* ));
557 
558  return ret;
559 }
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
#define MAXIT
Definition: mpr_numeric.cc:274
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
bool isfloat(gmp_complex **a)
Definition: mpr_numeric.cc:638
int k
Definition: cfEzgcd.cc:93
gmp_complex numbers based on
Definition: mpr_complex.h:178
#define WarnS
Definition: emacs.cc:81
#define omAlloc(size)
Definition: omAllocDecl.h:210
int j
Definition: myNF.cc:70
void divquad(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:671
void sortroots(gmp_complex **roots, int r, int c, bool isf)
Definition: mpr_numeric.cc:748
#define ST_ROOTS_LGSTEP
Definition: mpr_global.h:80
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:563
int i
Definition: cfEzgcd.cc:123
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
Definition: mpr_numeric.cc:695
Variable x
Definition: cfModGcd.cc:4023
void divlin(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:651
gmp_complex& rootContainer::operator[] ( const int  i)
inline

Definition at line 82 of file mpr_numeric.h.

83  {
84  return *theroots[i];
85  }
int i
Definition: cfEzgcd.cc:123
gmp_complex ** theroots
Definition: mpr_numeric.h:139
void rootContainer::solvequad ( gmp_complex **  a,
gmp_complex **  r,
int &  k,
int &  j 
)
private

Definition at line 695 of file mpr_numeric.cc.

696 {
697  gmp_float zero(0.0);
698 
699  if ((j>k)
700  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
701  {
702  gmp_complex sq(zero);
703  gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
704  gmp_complex disk((h1 * h1) - h2);
705  if (disk.imag().isZero())
706  {
707  if (disk.real()<zero)
708  {
709  sq.real(zero);
710  sq.imag(sqrt(-disk.real()));
711  }
712  else
713  sq = (gmp_complex)sqrt(disk.real());
714  }
715  else
716  sq = sqrt(disk);
717  *r[k+1] = sq - h1;
718  sq += h1;
719  *r[k] = (gmp_complex)0.0-sq;
720  if(sq.imag().isZero())
721  {
722  k = j;
723  j++;
724  }
725  else
726  {
727  j = k;
728  k--;
729  }
730  }
731  else
732  {
733  if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero()))
734  {
735  WerrorS("precision lost, try again with higher precision");
736  }
737  else
738  {
739  *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
740  if(r[k]->imag().isZero())
741  j++;
742  else
743  k--;
744  }
745  }
746 }
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
gmp_complex numbers based on
Definition: mpr_complex.h:178
int j
Definition: myNF.cc:70
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:329
bool isZero()
Definition: mpr_complex.h:241
bool isZero(const CFArray &A)
checks if entries of A are zero
bool rootContainer::solver ( const int  polishmode = PM_NONE)

Definition at line 450 of file mpr_numeric.cc.

451 {
452  int i;
453 
454  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
455  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
456  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
457 
458  // copy the coefficients of type number to type gmp_complex
459  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
460  for ( i=0; i <= tdg; i++ )
461  {
462  ad[i]= new gmp_complex();
463  if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i], currRing->cf );
464  }
465 
466  // now solve
467  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
468  if (!found_roots)
469  WarnS("rootContainer::solver: No roots found!");
470 
471  // free memory
472  for ( i=0; i <= tdg; i++ ) delete ad[i];
473  omFreeSize( (void *) ad, (tdg+1)*sizeof(gmp_complex*) );
474 
475  return found_roots;
476 }
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
gmp_complex numbers based on
Definition: mpr_complex.h:178
#define WarnS
Definition: emacs.cc:81
#define omAlloc(size)
Definition: omAllocDecl.h:210
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
gmp_complex numberToComplex(number num, const coeffs r)
Definition: mpr_complex.h:312
The main handler for Singular numbers which are suitable for Singular polynomials.
int i
Definition: cfEzgcd.cc:123
gmp_complex ** theroots
Definition: mpr_numeric.h:139
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:480
void rootContainer::sortre ( gmp_complex **  r,
int  l,
int  u,
int  inc 
)
private

Definition at line 767 of file mpr_numeric.cc.

768 {
769  int pos,i;
770  gmp_complex *x,*y;
771 
772  pos = l;
773  x = r[pos];
774  for (i=l+inc; i<=u; i+=inc)
775  {
776  if (r[i]->real()<x->real())
777  {
778  pos = i;
779  x = r[pos];
780  }
781  }
782  if (pos>l)
783  {
784  if (inc==1)
785  {
786  for (i=pos; i>l; i--)
787  r[i] = r[i-1];
788  r[l] = x;
789  }
790  else
791  {
792  y = r[pos+1];
793  for (i=pos+1; i+1>l; i--)
794  r[i] = r[i-2];
795  if (x->imag()>y->imag())
796  {
797  r[l] = x;
798  r[l+1] = y;
799  }
800  else
801  {
802  r[l] = y;
803  r[l+1] = x;
804  }
805  }
806  }
807  else if ((inc==2)&&(x->imag()<r[l+1]->imag()))
808  {
809  r[l] = r[l+1];
810  r[l+1] = x;
811  }
812 }
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
gmp_complex numbers based on
Definition: mpr_complex.h:178
int i
Definition: cfEzgcd.cc:123
gmp_float imag() const
Definition: mpr_complex.h:235
Variable x
Definition: cfModGcd.cc:4023
int l
Definition: cfEzgcd.cc:94
gmp_float real() const
Definition: mpr_complex.h:234
void rootContainer::sortroots ( gmp_complex **  roots,
int  r,
int  c,
bool  isf 
)
private

Definition at line 748 of file mpr_numeric.cc.

749 {
750  int j;
751 
752  for (j=0; j<r; j++) // the real roots
753  sortre(ro, j, r, 1);
754  if (c>=tdg) return;
755  if (isf)
756  {
757  for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly
758  sortre(ro, j, tdg-1, 2);
759  }
760  else
761  {
762  for (j=c; j+1<tdg; j++) // the complex roots for a general poly
763  sortre(ro, j, tdg-1, 1);
764  }
765 }
void sortre(gmp_complex **r, int l, int u, int inc)
Definition: mpr_numeric.cc:767
const ring r
Definition: syzextra.cc:208
int j
Definition: myNF.cc:70
bool rootContainer::swapRoots ( const int  from,
const int  to 
)

Definition at line 430 of file mpr_numeric.cc.

431 {
432  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
433  {
434  if ( to != from )
435  {
436  gmp_complex tmp( *theroots[from] );
437  *theroots[from]= *theroots[to];
438  *theroots[to]= tmp;
439  }
440  return true;
441  }
442 
443  // warning
444  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
445  return false;
446 }
gmp_complex numbers based on
Definition: mpr_complex.h:178
clock_t to
Definition: walk.cc:100
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define Warn
Definition: emacs.cc:80

Field Documentation

int rootContainer::anz
private

Definition at line 141 of file mpr_numeric.h.

number* rootContainer::coeffs
private

Definition at line 135 of file mpr_numeric.h.

bool rootContainer::found_roots
private

Definition at line 142 of file mpr_numeric.h.

number* rootContainer::ievpoint
private

Definition at line 136 of file mpr_numeric.h.

rootType rootContainer::rt
private

Definition at line 137 of file mpr_numeric.h.

int rootContainer::tdg
private

Definition at line 133 of file mpr_numeric.h.

gmp_complex** rootContainer::theroots
private

Definition at line 139 of file mpr_numeric.h.

int rootContainer::var
private

Definition at line 132 of file mpr_numeric.h.


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