 |
My Project
debian-1:4.1.2-p1+ds-2
|
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 64 of file mpr_numeric.h.
◆ rootType
Enumerator |
---|
none | |
cspecial | |
cspecialmu | |
det | |
onepoly | |
Definition at line 67 of file mpr_numeric.h.
◆ rootContainer() [1/2]
rootContainer::rootContainer |
( |
| ) |
|
◆ ~rootContainer()
rootContainer::~rootContainer |
( |
| ) |
|
◆ rootContainer() [2/2]
◆ checkimag()
◆ computefx()
Definition at line 799 of file mpr_numeric.cc.
812 for (
k=
m-1;
k >= 0;
k-- )
814 f2 = (
x * f2 ) + f1;
815 f1 = (
x * f1 ) + f0;
816 f0 = (
x * f0 ) + *a[
k];
817 ef =
abs( f0 ) + ( ex * ef );
◆ computegx()
Definition at line 820 of file mpr_numeric.cc.
833 for (
k= 1;
k <=
m;
k++ )
835 f2 = (
x * f2 ) + f1;
836 f1 = (
x * f1 ) + f0;
837 f0 = (
x * f0 ) + *a[
k];
838 ef =
abs( f0 ) + ( ex * ef );
◆ divlin()
Definition at line 636 of file mpr_numeric.cc.
644 for (
i=
j-1;
i > 0;
i-- )
645 *a[
i] += (*a[
i+1]*
x);
646 for (
i= 0;
i <
j;
i++ )
652 for (
i= 1;
i <
j;
i++)
653 *a[
i] += (*a[
i-1]*
y);
◆ divquad()
Definition at line 656 of file mpr_numeric.cc.
661 q((
x.real()*
x.real())+(
x.imag()*
x.imag()));
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++ )
676 for (
i= 2;
i <
j-1;
i++)
677 *a[
i] += ((*a[
i-1]*
p)-(*a[
i-2]*q));
◆ evPointCoord()
Definition at line 386 of file mpr_numeric.cc.
389 if (! ((
i >= 0) && (
i <
anz+2) ) )
390 WarnS(
"rootContainer::evPointCoord: index out of range");
392 WarnS(
"rootContainer::evPointCoord: ievpoint == NULL");
404 Warn(
"rootContainer::evPointCoord: NULL index %d",
i);
409 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 548 of file mpr_numeric.cc.
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 };
573 if ((fabs==zero) || (
abs(d)==zero))
return;
580 h= g2 - (((
f+
f) /
b ));
581 sq=
sqrt(( (
h * deg ) - g2 ) * (deg - one));
590 if((
gp.real()==zero)&&(
gp.imag()==zero))
603 if (*
x == x1)
goto ende;
607 if (
j %
MT ) *
x= x1;
608 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 465 of file mpr_numeric.cc.
471 bool ret=
true, isf=
isfloat(a), type=
true;
494 WarnS(
"Laguerre solver: Too many iterations!");
503 WarnS(
"Laguerre solver: Too many iterations in polish!");
508 if ((!type)&&(!((
x.real()==zero)&&(
x.imag()==zero))))
x = o/
x;
509 if (
x.imag() == zero)
541 for (
i=0;
i <=
tdg;
i++ )
delete ad[
i];
◆ operator[]()
◆ solvequad()
Definition at line 680 of file mpr_numeric.cc.
686 &&((!(*a[2]).real().
isZero())||(!(*a[2]).imag().
isZero())))
689 gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
691 if (disk.imag().isZero())
693 if (disk.real()<zero)
696 sq.imag(
sqrt(-disk.real()));
706 if(sq.imag().isZero())
719 if (((*a[1]).real().
isZero()) && ((*a[1]).imag().isZero()))
721 WerrorS(
"precision lost, try again with higher precision");
726 if(r[
k]->imag().isZero())
◆ solver()
bool rootContainer::solver |
( |
const int |
polishmode = PM_NONE | ) |
|
Definition at line 435 of file mpr_numeric.cc.
446 for (
i=0;
i <=
tdg;
i++ )
455 WarnS(
"rootContainer::solver: No roots found!");
458 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 752 of file mpr_numeric.cc.
760 for (
i=
l+inc;
i<=u;
i+=inc)
762 if (r[
i]->real()<
x->real())
772 for (
i=pos;
i>
l;
i--)
779 for (
i=pos+1;
i+1>
l;
i--)
781 if (
x->imag()>
y->imag())
793 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 415 of file mpr_numeric.cc.
418 if (
found_roots && ( from >= 0) && ( from <
tdg ) && ( to >= 0) && ( to <
tdg ) )
430 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:
EXTERN_VAR size_t gmp_output_digits
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)
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)
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
void sortroots(gmp_complex **roots, int r, int c, bool isf)
gmp_float cos(const gmp_float &a)
gmp_complex numbers based on