Loading...
Searching...
No Matches
mpr_numeric.h
Go to the documentation of this file.
#ifndef MPR_NUMERIC_H
#define MPR_NUMERIC_H
/****************************************
* Computer Algebra System SINGULAR *
****************************************/
/*
* ABSTRACT - multipolynomial resultants - numeric stuff
* ( root finder, vandermonde system solver, simplex )
*/
//-> include & define stuff
#include "coeffs/numbers.h"
#include "coeffs/mpr_global.h"
#include "coeffs/mpr_complex.h"
// define polish mode when finding roots
#define PM_NONE 0
#define PM_POLISH 1
#define PM_CORRUPT 2
//<-
//-> vandermonde system solver
/**
* vandermonde system solver for interpolating polynomials from their values
*/
class vandermonde
{
public:
vandermonde( const long _cn, const long _n,
const long _maxdeg, number *_p, const bool _homog = true );
~vandermonde();
/** Solves the Vandermode linear system
* \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
* Any computations are done using type number to get high pecision results.
* @param q n-tuple of results (right hand of equations)
* @return w n-tuple of coefficients of resulting polynomial, lowest deg first
*/
number * interpolateDense( const number * q );
poly numvec2poly(const number * q );
private:
void init();
private:
long n; // number of variables
long cn; // real number of coefficients of poly to interpolate
long maxdeg; // degree of the polynomial to interpolate
long l; // max number of coefficients in poly of deg maxdeg = (maxdeg+1)^n
number *p; // evaluation point
number *x; // coefficients, determinend by init() from *p
bool homog;
};
//<-
//-> rootContainer
/**
* complex root finder for univariate polynomials based on laguers algorithm
*/
class rootContainer
{
public:
enum rootType { none, cspecial, cspecialmu, det, onepoly };
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_complex & operator[] ( const int i );
inline gmp_complex & operator[] ( const int i )
{
return *theroots[i];
}
gmp_complex & evPointCoord( const int i );
inline gmp_complex * getRoot( const int i )
{
return theroots[i];
}
bool swapRoots( const int from, const int to );
int getAnzElems() { return anz; }
int getLDim() { return anz; }
int getAnzRoots() { return tdg; }
private:
rootContainer( const rootContainer & v );
/** Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg]
* (generated from the number coefficients coeffs[0..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.
*/
bool laguer_driver( gmp_complex ** a, gmp_complex ** roots, bool polish = true );
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);
/** 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.
*/
void laguer(gmp_complex ** a, int m, gmp_complex * x, int * its, bool type);
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);
int var;
int tdg;
number * coeffs;
number * ievpoint;
rootType rt;
gmp_complex ** theroots;
int anz;
bool found_roots;
};
//<-
class slists; typedef slists * lists;
//-> class rootArranger
class rootArranger
{
public:
friend lists listOfRoots( rootArranger*, const unsigned int oprec );
rootArranger( rootContainer ** _roots,
rootContainer ** _mu,
const int _howclean = PM_CORRUPT );
~rootArranger() {}
void solve_all();
void arrange();
bool success() { return found_roots; }
private:
rootArranger( const rootArranger & );
rootContainer ** roots;
rootContainer ** mu;
int howclean;
int rc,mc;
bool found_roots;
};
//<-
//-> simplex computation
// (used by sparse matrix construction)
#define SIMPLEX_EPS 1.0e-12
/** Linear Programming / Linear Optimization using Simplex - Algorithm
*
* On output, the tableau LiPM is indexed by two arrays of integers.
* ipsov[j] contains, for j=1..m, the number i whose original variable
* is now represented by row j+1 of LiPM. (left-handed vars in solution)
* (first row is the one with the objective function)
* izrov[j] contains, for j=1..n, the number i whose original variable
* x_i is now a right-handed variable, rep. by column j+1 of LiPM.
* These vars are all zero in the solution. The meaning of n<i<n+m1+m2
* is the same as above.
*/
class simplex
{
public:
int m; // number of constraints, make sure m == m1 + m2 + m3 !!
int n; // # of independent variables
int m1,m2,m3; // constraints <=, >= and ==
int icase; // == 0: finite solution found;
// == +1 objective function unbound; == -1: no solution
int *izrov,*iposv;
mprfloat **LiPM; // the matrix (of size [m+2, n+1])
/** #rows should be >= m+2, #cols >= n+1
*/
simplex( int rows, int cols );
~simplex();
BOOLEAN mapFromMatrix( matrix m );
matrix mapToMatrix( matrix m );
intvec * posvToIV();
intvec * zrovToIV();
void compute();
private:
simplex( const simplex & );
void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax );
void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 );
void simp3( mprfloat **a, int i1, int k1, int ip, int kp );
int LiPM_cols,LiPM_rows;
};
//<-
#endif /*MPR_NUMERIC_H*/
// local Variables: ***
// folded-file: t ***
// compile-command-1: "make installg" ***
// compile-command-2: "make install" ***
// End: ***
BOOLEAN
int BOOLEAN
Definition auxiliary.h:88
l
int l
Definition cfEzgcd.cc:100
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
x
Variable x
Definition cfModGcd.cc:4090
gmp_complex
gmp_complex numbers based on
Definition mpr_complex.h:181
gmp_float
Definition mpr_complex.h:32
intvec
Definition intvec.h:24
rootArranger::rc
int rc
Definition mpr_numeric.h:171
rootArranger::solve_all
void solve_all()
Definition mpr_numeric.cc:858
rootArranger::roots
rootContainer ** roots
Definition mpr_numeric.h:167
rootArranger::listOfRoots
friend lists listOfRoots(rootArranger *, const unsigned int oprec)
Definition ipshell.cc:5082
rootArranger::rootArranger
rootArranger(rootContainer **_roots, rootContainer **_mu, const int _howclean=PM_CORRUPT)
Definition mpr_numeric.cc:848
rootArranger::rootArranger
rootArranger(const rootArranger &)
rootArranger::mc
int mc
Definition mpr_numeric.h:171
rootArranger::~rootArranger
~rootArranger()
Definition mpr_numeric.h:157
rootArranger::mu
rootContainer ** mu
Definition mpr_numeric.h:168
rootArranger::found_roots
bool found_roots
Definition mpr_numeric.h:172
rootArranger::howclean
int howclean
Definition mpr_numeric.h:170
rootArranger::success
bool success()
Definition mpr_numeric.h:162
rootArranger::arrange
void arrange()
Definition mpr_numeric.cc:883
rootContainer
complex root finder for univariate polynomials based on laguers algorithm
Definition mpr_numeric.h:66
rootContainer::sortre
void sortre(gmp_complex **r, int l, int u, int inc)
Definition mpr_numeric.cc:754
rootContainer::getPoly
poly getPoly()
Definition mpr_numeric.cc:334
rootContainer::operator[]
gmp_complex & operator[](const int i)
Definition mpr_numeric.h:82
rootContainer::~rootContainer
~rootContainer()
Definition mpr_numeric.cc:277
rootContainer::laguer_driver
bool laguer_driver(gmp_complex **a, gmp_complex **roots, bool polish=true)
Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg] (generated from the number coeffic...
Definition mpr_numeric.cc:467
rootContainer::computegx
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:822
rootContainer::getRoot
gmp_complex * getRoot(const int i)
Definition mpr_numeric.h:88
rootContainer::fillContainer
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition mpr_numeric.cc:300
rootContainer::laguer
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:550
rootContainer::divlin
void divlin(gmp_complex **a, gmp_complex x, int j)
Definition mpr_numeric.cc:638
rootContainer::sortroots
void sortroots(gmp_complex **roots, int r, int c, bool isf)
Definition mpr_numeric.cc:735
rootContainer::rt
rootType rt
Definition mpr_numeric.h:137
rootContainer::divquad
void divquad(gmp_complex **a, gmp_complex x, int j)
Definition mpr_numeric.cc:658
rootContainer::swapRoots
bool swapRoots(const int from, const int to)
Definition mpr_numeric.cc:417
rootContainer::theroots
gmp_complex ** theroots
Definition mpr_numeric.h:139
rootContainer::computefx
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:801
rootContainer::coeffs
number * coeffs
Definition mpr_numeric.h:135
rootContainer::rootContainer
rootContainer()
Definition mpr_numeric.cc:264
rootContainer::solvequad
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
Definition mpr_numeric.cc:682
rootContainer::rootType
rootType
Definition mpr_numeric.h:68
rootContainer::onepoly
@ onepoly
Definition mpr_numeric.h:68
rootContainer::cspecial
@ cspecial
Definition mpr_numeric.h:68
rootContainer::none
@ none
Definition mpr_numeric.h:68
rootContainer::det
@ det
Definition mpr_numeric.h:68
rootContainer::cspecialmu
@ cspecialmu
Definition mpr_numeric.h:68
rootContainer::evPointCoord
gmp_complex & evPointCoord(const int i)
Definition mpr_numeric.cc:388
rootContainer::found_roots
bool found_roots
Definition mpr_numeric.h:142
rootContainer::isfloat
bool isfloat(gmp_complex **a)
Definition mpr_numeric.cc:625
rootContainer::rootContainer
rootContainer(const rootContainer &v)
rootContainer::getLDim
int getLDim()
Definition mpr_numeric.h:96
rootContainer::checkimag
void checkimag(gmp_complex *x, gmp_float &e)
Definition mpr_numeric.cc:617
rootContainer::anz
int anz
Definition mpr_numeric.h:141
rootContainer::getAnzRoots
int getAnzRoots()
Definition mpr_numeric.h:97
rootContainer::solver
bool solver(const int polishmode=PM_NONE)
Definition mpr_numeric.cc:437
rootContainer::ievpoint
number * ievpoint
Definition mpr_numeric.h:136
rootContainer::tdg
int tdg
Definition mpr_numeric.h:133
rootContainer::getAnzElems
int getAnzElems()
Definition mpr_numeric.h:95
rootContainer::var
int var
Definition mpr_numeric.h:132
simplex::zrovToIV
intvec * zrovToIV()
Definition mpr_numeric.cc:1084
simplex::LiPM
mprfloat ** LiPM
Definition mpr_numeric.h:205
simplex::m1
int m1
Definition mpr_numeric.h:200
simplex::n
int n
Definition mpr_numeric.h:199
simplex::m2
int m2
Definition mpr_numeric.h:200
simplex::iposv
int * iposv
Definition mpr_numeric.h:203
simplex::LiPM_rows
int LiPM_rows
Definition mpr_numeric.h:225
simplex::mapFromMatrix
BOOLEAN mapFromMatrix(matrix m)
Definition mpr_numeric.cc:1011
simplex::izrov
int * izrov
Definition mpr_numeric.h:203
simplex::~simplex
~simplex()
Definition mpr_numeric.cc:997
simplex::icase
int icase
Definition mpr_numeric.h:201
simplex::simplex
simplex(const simplex &)
simplex::compute
void compute()
Definition mpr_numeric.cc:1095
simplex::LiPM_cols
int LiPM_cols
Definition mpr_numeric.h:225
simplex::simplex
simplex(int rows, int cols)
rows should be >= m+2, cols >= n+1
Definition mpr_numeric.cc:972
simplex::m3
int m3
Definition mpr_numeric.h:200
simplex::simp2
void simp2(mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
Definition mpr_numeric.cc:1298
simplex::m
int m
Definition mpr_numeric.h:198
simplex::simp3
void simp3(mprfloat **a, int i1, int k1, int ip, int kp)
Definition mpr_numeric.cc:1337
simplex::simp1
void simp1(mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
Definition mpr_numeric.cc:1263
simplex::mapToMatrix
matrix mapToMatrix(matrix m)
Definition mpr_numeric.cc:1040
simplex::posvToIV
intvec * posvToIV()
Definition mpr_numeric.cc:1073
slists
Definition lists.h:24
vandermonde::init
void init()
Definition mpr_numeric.cc:53
vandermonde::homog
bool homog
Definition mpr_numeric.h:57
vandermonde::numvec2poly
poly numvec2poly(const number *q)
Definition mpr_numeric.cc:93
vandermonde::x
number * x
Definition mpr_numeric.h:55
vandermonde::vandermonde
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)
Definition mpr_numeric.cc:35
vandermonde::p
number * p
Definition mpr_numeric.h:54
vandermonde::maxdeg
long maxdeg
Definition mpr_numeric.h:51
vandermonde::l
long l
Definition mpr_numeric.h:52
vandermonde::n
long n
Definition mpr_numeric.h:49
vandermonde::interpolateDense
number * interpolateDense(const number *q)
Solves the Vandermode linear system \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
Definition mpr_numeric.cc:146
vandermonde::cn
long cn
Definition mpr_numeric.h:50
vandermonde::~vandermonde
~vandermonde()
Definition mpr_numeric.cc:46
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
matrix
ip_smatrix * matrix
Definition matpol.h:43
mpr_complex.h
mpr_global.h
mprfloat
double mprfloat
Definition mpr_global.h:17
PM_NONE
#define PM_NONE
Definition mpr_numeric.h:19
PM_CORRUPT
#define PM_CORRUPT
Definition mpr_numeric.h:21
lists
slists * lists
Definition mpr_numeric.h:146
numbers.h