 |
My Project
UNKNOWN_GIT_VERSION
|
complex root finder for univariate polynomials based on laguers algorithm
More...
#include <mpr_numeric.h>
|
| 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) |
|
complex root finder for univariate polynomials based on laguers algorithm
Definition at line 65 of file mpr_numeric.h.
◆ rootType
Enumerator |
---|
none | |
cspecial | |
cspecialmu | |
det | |
onepoly | |
Definition at line 68 of file mpr_numeric.h.
◆ rootContainer() [1/2]
rootContainer::rootContainer |
( |
| ) |
|
◆ ~rootContainer()
rootContainer::~rootContainer |
( |
| ) |
|
◆ rootContainer() [2/2]
◆ checkimag()
◆ computefx()
Definition at line 805 of file mpr_numeric.cc.
817 for (
k=
m-1;
k >= 0;
k-- )
819 f2 = (
x * f2 ) + f1;
820 f1 = (
x * f1 ) + f0;
821 f0 = (
x * f0 ) + *a[
k];
822 ef =
abs( f0 ) + ( ex * ef );
◆ computegx()
Definition at line 826 of file mpr_numeric.cc.
838 for (
k= 1;
k <=
m;
k++ )
840 f2 = (
x * f2 ) + f1;
841 f1 = (
x * f1 ) + f0;
842 f0 = (
x * f0 ) + *a[
k];
843 ef =
abs( f0 ) + ( ex * ef );
◆ divlin()
Definition at line 642 of file mpr_numeric.cc.
649 for (
i=
j-1;
i > 0;
i-- )
650 *a[
i] += (*a[
i+1]*
x);
651 for (
i= 0;
i <
j;
i++ )
657 for (
i= 1;
i <
j;
i++)
658 *a[
i] += (*a[
i-1]*
y);
◆ divquad()
Definition at line 662 of file mpr_numeric.cc.
666 q((
x.real()*
x.real())+(
x.imag()*
x.imag()));
670 *a[
j-1] += (*a[
j]*
p);
671 for (
i=
j-2;
i > 1;
i-- )
672 *a[
i] += ((*a[
i+1]*
p)-(*a[
i+2]*q));
673 for (
i= 0;
i <
j-1;
i++ )
681 for (
i= 2;
i <
j-1;
i++)
682 *a[
i] += ((*a[
i-1]*
p)-(*a[
i-2]*q));
◆ evPointCoord()
Definition at line 392 of file mpr_numeric.cc.
394 if (! ((
i >= 0) && (
i <
anz+2) ) )
395 WarnS(
"rootContainer::evPointCoord: index out of range");
397 WarnS(
"rootContainer::evPointCoord: ievpoint == NULL");
409 Warn(
"rootContainer::evPointCoord: NULL index %d",
i);
414 Warn(
"rootContainer::evPointCoord: Wrong index %d, found_roots %s",
i,
found_roots?
"true":
"false");
◆ fillContainer()
void rootContainer::fillContainer |
( |
number * |
_coeffs, |
|
|
number * |
_ievpoint, |
|
|
const int |
_var, |
|
|
const int |
_tdg, |
|
|
const rootType |
_rt, |
|
|
const int |
_anz |
|
) |
| |
◆ getAnzElems()
int rootContainer::getAnzElems |
( |
| ) |
|
|
inline |
◆ getAnzRoots()
int rootContainer::getAnzRoots |
( |
| ) |
|
|
inline |
◆ getLDim()
int rootContainer::getLDim |
( |
| ) |
|
|
inline |
◆ getPoly()
poly rootContainer::getPoly |
( |
| ) |
|
◆ getRoot()
◆ isfloat()
◆ laguer()
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 554 of file mpr_numeric.cc.
559 gmp_complex dx, x1,
b, d,
f,
g,
h, sq,
gp, gm, g2;
560 gmp_float frac_g[
MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
578 if ((fabs==zero) || (
abs(d)==zero))
return;
585 h= g2 - (((
f+
f) /
b ));
586 sq=
sqrt(( (
h * deg ) - g2 ) * (deg - one));
595 if((
gp.real()==zero)&&(
gp.imag()==zero))
608 if (*
x == x1)
goto ende;
612 if (
j %
MT ) *
x= x1;
613 else *
x -= ( dx * frac_g[
j /
MT ] );
◆ laguer_driver()
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 471 of file mpr_numeric.cc.
476 bool ret=
true, isf=
isfloat(a), type=
true;
499 WarnS(
"Laguerre solver: Too many iterations!");
508 WarnS(
"Laguerre solver: Too many iterations in polish!");
513 if ((!type)&&(!((
x.real()==zero)&&(
x.imag()==zero))))
x = o/
x;
514 if (
x.imag() == zero)
546 for (
i=0;
i <=
tdg;
i++ )
delete ad[
i];
◆ operator[]()
◆ solvequad()
Definition at line 686 of file mpr_numeric.cc.
691 &&((!(*a[2]).real().
isZero())||(!(*a[2]).imag().
isZero())))
694 gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
696 if (disk.imag().isZero())
698 if (disk.real()<zero)
701 sq.imag(
sqrt(-disk.real()));
711 if(sq.imag().isZero())
724 if (((*a[1]).real().
isZero()) && ((*a[1]).imag().isZero()))
726 WerrorS(
"precision lost, try again with higher precision");
731 if(r[
k]->imag().isZero())
◆ solver()
bool rootContainer::solver |
( |
const int |
polishmode = PM_NONE | ) |
|
Definition at line 441 of file mpr_numeric.cc.
451 for (
i=0;
i <=
tdg;
i++ )
460 WarnS(
"rootContainer::solver: No roots found!");
463 for (
i=0;
i <=
tdg;
i++ )
delete ad[
i];
◆ sortre()
void rootContainer::sortre |
( |
gmp_complex ** |
r, |
|
|
int |
l, |
|
|
int |
u, |
|
|
int |
inc |
|
) |
| |
|
private |
Definition at line 758 of file mpr_numeric.cc.
765 for (
i=
l+inc;
i<=u;
i+=inc)
767 if (r[
i]->real()<
x->real())
777 for (
i=pos;
i>
l;
i--)
784 for (
i=pos+1;
i+1>
l;
i--)
786 if (
x->imag()>
y->imag())
798 else if ((inc==2)&&(
x->imag()<r[
l+1]->imag()))
◆ sortroots()
void rootContainer::sortroots |
( |
gmp_complex ** |
roots, |
|
|
int |
r, |
|
|
int |
c, |
|
|
bool |
isf |
|
) |
| |
|
private |
◆ swapRoots()
bool rootContainer::swapRoots |
( |
const int |
from, |
|
|
const int |
to |
|
) |
| |
Definition at line 421 of file mpr_numeric.cc.
423 if (
found_roots && ( from >= 0) && ( from <
tdg ) && ( to >= 0) && ( to <
tdg ) )
435 Warn(
" rootContainer::changeRoots: Wrong index %d, %d",from,to);
◆ anz
◆ coeffs
number* rootContainer::coeffs |
|
private |
◆ found_roots
bool rootContainer::found_roots |
|
private |
◆ ievpoint
number* rootContainer::ievpoint |
|
private |
◆ rt
◆ tdg
◆ theroots
◆ var
The documentation for this class was generated from the following files:
bool isZero(const CFArray &A)
checks if entries of A are zero
void divlin(gmp_complex **a, gmp_complex x, int j)
const CanonicalForm int const CFList const Variable & y
bool isfloat(gmp_complex **a)
gmp_float sqrt(const gmp_float &a)
gmp_float sin(const gmp_float &a)
#define mprSTICKYPROT(msg)
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Rational abs(const Rational &a)
#define omFreeSize(addr, size)
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 ...
void checkimag(gmp_complex *x, gmp_float &e)
void sortre(gmp_complex **r, int l, int u, int inc)
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...
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
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)
gmp_complex numberToComplex(number num, const coeffs r)
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
void divquad(gmp_complex **a, gmp_complex x, int j)
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 WerrorS(const char *s)
void sortroots(gmp_complex **roots, int r, int c, bool isf)
gmp_float cos(const gmp_float &a)
gmp_complex numbers based on