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 277 of file mpr_numeric.cc.

278 {
279  rt=none;
280 
281  coeffs= NULL;
282  ievpoint= NULL;
283  theroots= NULL;
284 
285  found_roots= false;
286 }
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::~rootContainer ( )

Definition at line 290 of file mpr_numeric.cc.

291 {
292  int i;
293  // free coeffs, ievpoint
294  if ( ievpoint != NULL )
295  {
296  for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
297  omFreeSize( (void *)ievpoint, (anz+2) * sizeof( number ) );
298  }
299 
300  for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i );
301  omFreeSize( (void *)coeffs, (tdg+1) * sizeof( number ) );
302 
303  // theroots löschen
304  for ( i=0; i < tdg; i++ ) delete theroots[i];
305  omFreeSize( (void *) theroots, (tdg)*sizeof(gmp_complex*) );
306 
307  //mprPROTnl("~rootContainer()");
308 }
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() [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 629 of file mpr_numeric.cc.

630 {
631  if(abs(x->imag())<abs(x->real())*e)
632  {
633  x->imag(0.0);
634  }
635 }
gmp_float real() const
Definition: mpr_complex.h:234
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
gmp_float imag() const
Definition: mpr_complex.h:235

§ 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 813 of file mpr_numeric.cc.

816 {
817  int k;
818 
819  f0= *a[m];
820  ef= abs(f0);
821  f1= gmp_complex(0.0);
822  f2= f1;
823  ex= abs(x);
824 
825  for ( k= m-1; k >= 0; k-- )
826  {
827  f2 = ( x * f2 ) + f1;
828  f1 = ( x * f1 ) + f0;
829  f0 = ( x * f0 ) + *a[k];
830  ef = abs( f0 ) + ( ex * ef );
831  }
832 }
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

§ 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 834 of file mpr_numeric.cc.

837 {
838  int k;
839 
840  f0= *a[0];
841  ef= abs(f0);
842  f1= gmp_complex(0.0);
843  f2= f1;
844  ex= abs(x);
845 
846  for ( k= 1; k <= m; k++ )
847  {
848  f2 = ( x * f2 ) + f1;
849  f1 = ( x * f1 ) + f0;
850  f0 = ( x * f0 ) + *a[k];
851  ef = abs( f0 ) + ( ex * ef );
852  }
853 }
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

§ divlin()

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

Definition at line 650 of file mpr_numeric.cc.

651 {
652  int i;
653  gmp_float o(1.0);
654 
655  if (abs(x)<o)
656  {
657  for (i= j-1; i > 0; i-- )
658  *a[i] += (*a[i+1]*x);
659  for (i= 0; i < j; i++ )
660  *a[i] = *a[i+1];
661  }
662  else
663  {
664  gmp_complex y(o/x);
665  for (i= 1; i < j; i++)
666  *a[i] += (*a[i-1]*y);
667  }
668 }
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

§ divquad()

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

Definition at line 670 of file mpr_numeric.cc.

671 {
672  int i;
673  gmp_float o(1.0),p(x.real()+x.real()),
674  q((x.real()*x.real())+(x.imag()*x.imag()));
675 
676  if (abs(x)<o)
677  {
678  *a[j-1] += (*a[j]*p);
679  for (i= j-2; i > 1; i-- )
680  *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
681  for (i= 0; i < j-1; i++ )
682  *a[i] = *a[i+2];
683  }
684  else
685  {
686  p = p/q;
687  q = o/q;
688  *a[1] += (*a[0]*p);
689  for (i= 2; i < j-1; i++)
690  *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
691  }
692 }
return P p
Definition: myNF.cc:203
gmp_float real() const
Definition: mpr_complex.h:234
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

§ evPointCoord()

gmp_complex & rootContainer::evPointCoord ( const int  i)

Definition at line 400 of file mpr_numeric.cc.

401 {
402  if (! ((i >= 0) && (i < anz+2) ) )
403  WarnS("rootContainer::evPointCoord: index out of range");
404  if (ievpoint == NULL)
405  WarnS("rootContainer::evPointCoord: ievpoint == NULL");
406 
407  if ( (rt == cspecialmu) && found_roots ) // FIX ME
408  {
409  if ( ievpoint[i] != NULL )
410  {
411  gmp_complex *tmp= new gmp_complex();
412  *tmp= numberToComplex(ievpoint[i], currRing->cf);
413  return *tmp;
414  }
415  else
416  {
417  Warn("rootContainer::evPointCoord: NULL index %d",i);
418  }
419  }
420 
421  // warning
422  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
423  gmp_complex *tmp= new gmp_complex();
424  return *tmp;
425 }
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:10
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

§ fillContainer()

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

Definition at line 312 of file mpr_numeric.cc.

315 {
316  int i;
317  number nn= nInit(0);
318  var=_var;
319  tdg=_tdg;
320  coeffs=_coeffs;
321  rt=_rt;
322  anz=_anz;
323 
324  for ( i=0; i <= tdg; i++ )
325  {
326  if ( nEqual(coeffs[i],nn) )
327  {
328  nDelete( &coeffs[i] );
329  coeffs[i]=NULL;
330  }
331  }
332  nDelete( &nn );
333 
334  if ( rt == cspecialmu && _ievpoint ) // copy ievpoint
335  {
336  ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
337  for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
338  }
339 
340  theroots= NULL;
341  found_roots= false;
342 }
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

§ 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 346 of file mpr_numeric.cc.

347 {
348  int i;
349 
350  poly result= NULL;
351  poly ppos;
352 
353  if ( (rt == cspecial) || ( rt == cspecialmu ) )
354  {
355  for ( i= tdg; i >= 0; i-- )
356  {
357  if ( coeffs[i] )
358  {
359  poly p= pOne();
360  //pSetExp( p, var+1, i);
361  pSetExp( p, 1, i);
362  pSetCoeff( p, nCopy( coeffs[i] ) );
363  pSetm( p );
364  if (result)
365  {
366  ppos->next=p;
367  ppos=ppos->next;
368  }
369  else
370  {
371  result=p;
372  ppos=p;
373  }
374 
375  }
376  }
377  if (result!=NULL) pSetm( result );
378  }
379 
380  return result;
381 }
#define pSetm(p)
Definition: polys.h:253
#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:297
#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

§ getRoot()

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

§ isfloat()

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

Definition at line 637 of file mpr_numeric.cc.

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

§ 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 562 of file mpr_numeric.cc.

563 {
564  int iter,j;
565  gmp_float zero(0.0),one(1.0),deg(m);
566  gmp_float abx_g, err_g, fabs;
567  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
568  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
569 
570  gmp_float epss(0.1);
571  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),gmp_output_digits);
572 
573  for ( iter= 1; iter <= MAXIT; iter++ )
574  {
576  *its=iter;
577  if (type)
578  computefx(a,*x,m,b,d,f,abx_g,err_g);
579  else
580  computegx(a,*x,m,b,d,f,abx_g,err_g);
581  err_g *= epss; // EPSS;
582 
583  fabs = abs(b);
584  if (fabs <= err_g)
585  {
586  if ((fabs==zero) || (abs(d)==zero)) return;
587  *x -= (b/d); // a last newton-step
588  goto ende;
589  }
590 
591  g= d / b;
592  g2 = g * g;
593  h= g2 - (((f+f) / b ));
594  sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
595  gp= g + sq;
596  gm= g - sq;
597  if (abs(gp)<abs(gm))
598  {
599  dx = deg/gm;
600  }
601  else
602  {
603  if((gp.real()==zero)&&(gp.imag()==zero))
604  {
605  dx.real(cos((mprfloat)iter));
606  dx.imag(sin((mprfloat)iter));
607  dx = dx*(one+abx_g);
608  }
609  else
610  {
611  dx = deg/gp;
612  }
613  }
614  x1= *x - dx;
615 
616  if (*x == x1) goto ende;
617 
618  j = iter%MMOD;
619  if (j==0) j=MT;
620  if ( j % MT ) *x= x1;
621  else *x -= ( dx * frac_g[ j / MT ] );
622  }
623 
624  *its= MAXIT+1;
625 ende:
626  checkimag(x,epss);
627 }
#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:834
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:813
#define MT
Definition: mpr_numeric.cc:271
f
Definition: cfModGcd.cc:4022
#define MAXIT
Definition: mpr_numeric.cc:273
CFFListIterator iter
Definition: facAbsBiFact.cc:54
void checkimag(gmp_complex *x, gmp_float &e)
Definition: mpr_numeric.cc:629
gmp_float real() const
Definition: mpr_complex.h:234
#define MMOD
Definition: mpr_numeric.cc:272
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
#define MR
Definition: mpr_numeric.cc:270
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 imag() const
Definition: mpr_complex.h:235

§ 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 479 of file mpr_numeric.cc.

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

§ operator[]()

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

§ solvequad()

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

Definition at line 694 of file mpr_numeric.cc.

695 {
696  gmp_float zero(0.0);
697 
698  if ((j>k)
699  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
700  {
701  gmp_complex sq(zero);
702  gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
703  gmp_complex disk((h1 * h1) - h2);
704  if (disk.imag().isZero())
705  {
706  if (disk.real()<zero)
707  {
708  sq.real(zero);
709  sq.imag(sqrt(-disk.real()));
710  }
711  else
712  sq = (gmp_complex)sqrt(disk.real());
713  }
714  else
715  sq = sqrt(disk);
716  *r[k+1] = sq - h1;
717  sq += h1;
718  *r[k] = (gmp_complex)0.0-sq;
719  if(sq.imag().isZero())
720  {
721  k = j;
722  j++;
723  }
724  else
725  {
726  j = k;
727  k--;
728  }
729  }
730  else
731  {
732  if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero()))
733  {
734  WerrorS("precision lost, try again with higher precision");
735  }
736  else
737  {
738  *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
739  if(r[k]->imag().isZero())
740  j++;
741  else
742  k--;
743  }
744  }
745 }
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

§ solver()

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

Definition at line 449 of file mpr_numeric.cc.

450 {
451  int i;
452 
453  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
454  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
455  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
456 
457  // copy the coefficients of type number to type gmp_complex
458  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
459  for ( i=0; i <= tdg; i++ )
460  {
461  ad[i]= new gmp_complex();
462  if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i], currRing->cf );
463  }
464 
465  // now solve
466  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
467  if (!found_roots)
468  WarnS("rootContainer::solver: No roots found!");
469 
470  // free memory
471  for ( i=0; i <= tdg; i++ ) delete ad[i];
472  omFreeSize( (void *) ad, (tdg+1)*sizeof(gmp_complex*) );
473 
474  return found_roots;
475 }
#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:10
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:479

§ sortre()

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

Definition at line 766 of file mpr_numeric.cc.

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

§ sortroots()

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

Definition at line 747 of file mpr_numeric.cc.

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

§ swapRoots()

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

Definition at line 429 of file mpr_numeric.cc.

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

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: