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 }
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: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++ ) nDelete( coeffs + i );
288  omFreeSize( (void *)coeffs, (tdg+1) * sizeof( number ) );
289 
290  // theroots löschen
291  for ( i=0; i < tdg; i++ ) delete theroots[i];
292  omFreeSize( (void *) theroots, (tdg)*sizeof(gmp_complex*) );
293 
294  //mprPROTnl("~rootContainer()");
295 }
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:125
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define nDelete(n)
Definition: numbers.h:16
#define NULL
Definition: omList.c:12

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

617 {
618  if(abs(x->imag())<abs(x->real())*e)
619  {
620  x->imag(0.0);
621  }
622 }
gmp_float real() const
Definition: mpr_complex.h:234
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
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 800 of file mpr_numeric.cc.

803 {
804  int k;
805 
806  f0= *a[m];
807  ef= abs(f0);
808  f1= gmp_complex(0.0);
809  f2= f1;
810  ex= abs(x);
811 
812  for ( k= m-1; k >= 0; k-- )
813  {
814  f2 = ( x * f2 ) + f1;
815  f1 = ( x * f1 ) + f0;
816  f0 = ( x * f0 ) + *a[k];
817  ef = abs( f0 ) + ( ex * ef );
818  }
819 }
int k
Definition: cfEzgcd.cc:92
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
int m
Definition: cfEzgcd.cc:121

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

824 {
825  int k;
826 
827  f0= *a[0];
828  ef= abs(f0);
829  f1= gmp_complex(0.0);
830  f2= f1;
831  ex= abs(x);
832 
833  for ( k= 1; k <= m; k++ )
834  {
835  f2 = ( x * f2 ) + f1;
836  f1 = ( x * f1 ) + f0;
837  f0 = ( x * f0 ) + *a[k];
838  ef = abs( f0 ) + ( ex * ef );
839  }
840 }
int k
Definition: cfEzgcd.cc:92
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
int m
Definition: cfEzgcd.cc:121

◆ divlin()

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

Definition at line 637 of file mpr_numeric.cc.

638 {
639  int i;
640  gmp_float o(1.0);
641 
642  if (abs(x)<o)
643  {
644  for (i= j-1; i > 0; i-- )
645  *a[i] += (*a[i+1]*x);
646  for (i= 0; i < j; i++ )
647  *a[i] = *a[i+1];
648  }
649  else
650  {
651  gmp_complex y(o/x);
652  for (i= 1; i < j; i++)
653  *a[i] += (*a[i-1]*y);
654  }
655 }
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
int j
Definition: facHensel.cc:105
gmp_complex numbers based on
Definition: mpr_complex.h:178
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
int i
Definition: cfEzgcd.cc:125

◆ divquad()

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

Definition at line 657 of file mpr_numeric.cc.

658 {
659  int i;
660  gmp_float o(1.0),p(x.real()+x.real()),
661  q((x.real()*x.real())+(x.imag()*x.imag()));
662 
663  if (abs(x)<o)
664  {
665  *a[j-1] += (*a[j]*p);
666  for (i= j-2; i > 1; i-- )
667  *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
668  for (i= 0; i < j-1; i++ )
669  *a[i] = *a[i+2];
670  }
671  else
672  {
673  p = p/q;
674  q = o/q;
675  *a[1] += (*a[0]*p);
676  for (i= 2; i < j-1; i++)
677  *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
678  }
679 }
int j
Definition: facHensel.cc:105
gmp_float real() const
Definition: mpr_complex.h:234
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
int i
Definition: cfEzgcd.cc:125
int p
Definition: cfModGcd.cc:4019
gmp_float imag() const
Definition: mpr_complex.h:235

◆ evPointCoord()

gmp_complex & rootContainer::evPointCoord ( const int  i)

Definition at line 387 of file mpr_numeric.cc.

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

◆ fillContainer()

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

Definition at line 299 of file mpr_numeric.cc.

302 {
303  int i;
304  number nn= nInit(0);
305  var=_var;
306  tdg=_tdg;
307  coeffs=_coeffs;
308  rt=_rt;
309  anz=_anz;
310 
311  for ( i=0; i <= tdg; i++ )
312  {
313  if ( nEqual(coeffs[i],nn) )
314  {
315  nDelete( &coeffs[i] );
316  coeffs[i]=NULL;
317  }
318  }
319  nDelete( &nn );
320 
321  if ( rt == cspecialmu && _ievpoint ) // copy ievpoint
322  {
323  ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
324  for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
325  }
326 
327  theroots= NULL;
328  found_roots= false;
329 }
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:125
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define nDelete(n)
Definition: numbers.h:16
#define NULL
Definition: omList.c:12
#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 333 of file mpr_numeric.cc.

334 {
335  int i;
336 
337  poly result= NULL;
338  poly ppos;
339 
340  if ( (rt == cspecial) || ( rt == cspecialmu ) )
341  {
342  for ( i= tdg; i >= 0; i-- )
343  {
344  if ( coeffs[i] )
345  {
346  poly p= pOne();
347  //pSetExp( p, var+1, i);
348  pSetExp( p, 1, i);
349  pSetCoeff( p, nCopy( coeffs[i] ) );
350  pSetm( p );
351  if (result)
352  {
353  ppos->next=p;
354  ppos=ppos->next;
355  }
356  else
357  {
358  result=p;
359  ppos=p;
360  }
361 
362  }
363  }
364  if (result!=NULL) pSetm( result );
365  }
366 
367  return result;
368 }
#define pSetm(p)
Definition: polys.h:266
#define pSetExp(p, i, v)
Definition: polys.h:42
rootType rt
Definition: mpr_numeric.h:137
The main handler for Singular numbers which are suitable for Singular polynomials.
int i
Definition: cfEzgcd.cc:125
#define pOne()
Definition: polys.h:310
#define NULL
Definition: omList.c:12
#define nCopy(n)
Definition: numbers.h:15
int p
Definition: cfModGcd.cc:4019
#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:125
gmp_complex ** theroots
Definition: mpr_numeric.h:139

◆ isfloat()

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

Definition at line 624 of file mpr_numeric.cc.

625 {
626  gmp_float z(0.0);
627  gmp_complex *b;
628  for (int i=tdg; i >= 0; i-- )
629  {
630  b = &(*a[i]);
631  if (!(b->imag()==z))
632  return false;
633  }
634  return true;
635 }
gmp_complex numbers based on
Definition: mpr_complex.h:178
CanonicalForm b
Definition: cfModGcd.cc:4044
int i
Definition: cfEzgcd.cc:125
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 549 of file mpr_numeric.cc.

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

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

◆ 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:125
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 681 of file mpr_numeric.cc.

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

437 {
438  int i;
439 
440  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
441  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
442  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
443 
444  // copy the coefficients of type number to type gmp_complex
445  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
446  for ( i=0; i <= tdg; i++ )
447  {
448  ad[i]= new gmp_complex();
449  if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i], currRing->cf );
450  }
451 
452  // now solve
453  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
454  if (!found_roots)
455  WarnS("rootContainer::solver: No roots found!");
456 
457  // free memory
458  for ( i=0; i <= tdg; i++ ) delete ad[i];
459  omFreeSize( (void *) ad, (tdg+1)*sizeof(gmp_complex*) );
460 
461  return found_roots;
462 }
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
gmp_complex numbers based on
Definition: mpr_complex.h:178
#define WarnS
Definition: emacs.cc:78
#define omAlloc(size)
Definition: omAllocDecl.h:210
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:125
gmp_complex ** theroots
Definition: mpr_numeric.h:139
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
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:466

◆ sortre()

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

Definition at line 753 of file mpr_numeric.cc.

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

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

◆ swapRoots()

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

Definition at line 416 of file mpr_numeric.cc.

417 {
418  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
419  {
420  if ( to != from )
421  {
422  gmp_complex tmp( *theroots[from] );
423  *theroots[from]= *theroots[to];
424  *theroots[to]= tmp;
425  }
426  return true;
427  }
428 
429  // warning
430  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
431  return false;
432 }
gmp_complex numbers based on
Definition: mpr_complex.h:178
gmp_complex ** theroots
Definition: mpr_numeric.h:139
#define Warn
Definition: emacs.cc:77

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: