My Project  debian-1:4.1.1-p2+ds-4build2
mpr_numeric.h
Go to the documentation of this file.
1 #ifndef MPR_NUMERIC_H
2 #define MPR_NUMERIC_H
3 /****************************************
4 * Computer Algebra System SINGULAR *
5 ****************************************/
6 
7 
8 /*
9 * ABSTRACT - multipolynomial resultants - numeric stuff
10 * ( root finder, vandermonde system solver, simplex )
11 */
12 
13 //-> include & define stuff
14 #include "coeffs/numbers.h"
15 #include "coeffs/mpr_global.h"
16 #include "coeffs/mpr_complex.h"
17 
18 // define polish mode when finding roots
19 #define PM_NONE 0
20 #define PM_POLISH 1
21 #define PM_CORRUPT 2
22 //<-
23 
24 //-> vandermonde system solver
25 /**
26  * vandermonde system solver for interpolating polynomials from their values
27  */
29 {
30 public:
31  vandermonde( const long _cn, const long _n,
32  const long _maxdeg, number *_p, const bool _homog = true );
33  ~vandermonde();
34 
35  /** Solves the Vandermode linear system
36  * \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
37  * Any computations are done using type number to get high pecision results.
38  * @param q n-tuple of results (right hand of equations)
39  * @return w n-tuple of coefficients of resulting polynomial, lowest deg first
40  */
41  number * interpolateDense( const number * q );
42 
43  poly numvec2poly(const number * q );
44 
45 private:
46  void init();
47 
48 private:
49  long n; // number of variables
50  long cn; // real number of coefficients of poly to interpolate
51  long maxdeg; // degree of the polynomial to interpolate
52  long l; // max number of coefficients in poly of deg maxdeg = (maxdeg+1)^n
53 
54  number *p; // evaluation point
55  number *x; // coefficients, determinend by init() from *p
56 
57  bool homog;
58 };
59 //<-
60 
61 //-> rootContainer
62 /**
63  * complex root finder for univariate polynomials based on laguers algorithm
64  */
66 {
67 public:
69 
70  rootContainer();
72 
73  void fillContainer( number *_coeffs, number *_ievpoint,
74  const int _var, const int _tdg,
75  const rootType _rt, const int _anz );
76 
77  bool solver( const int polishmode= PM_NONE );
78 
79  poly getPoly();
80 
81  //gmp_complex & operator[] ( const int i );
82  inline gmp_complex & operator[] ( const int i )
83  {
84  return *theroots[i];
85  }
86  gmp_complex & evPointCoord( const int i );
87 
88  inline gmp_complex * getRoot( const int i )
89  {
90  return theroots[i];
91  }
92 
93  bool swapRoots( const int from, const int to );
94 
95  int getAnzElems() { return anz; }
96  int getLDim() { return anz; }
97  int getAnzRoots() { return tdg; }
98 
99 private:
100  rootContainer( const rootContainer & v );
101 
102  /** Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg]
103  * (generated from the number coefficients coeffs[0..tdg]) of the polynomial
104  * this routine successively calls "laguer" and finds all m complex roots in
105  * roots[0..tdg]. The bool var "polish" should be input as "true" if polishing
106  * (also by "laguer") is desired, "false" if the roots will be subsequently
107  * polished by other means.
108  */
109  bool laguer_driver( gmp_complex ** a, gmp_complex ** roots, bool polish = true );
110  bool isfloat(gmp_complex **a);
111  void divlin(gmp_complex **a, gmp_complex x, int j);
112  void divquad(gmp_complex **a, gmp_complex x, int j);
113  void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j);
114  void sortroots(gmp_complex **roots, int r, int c, bool isf);
115  void sortre(gmp_complex **r, int l, int u, int inc);
116 
117  /** Given the degree m and the m+1 complex coefficients a[0..m] of the
118  * polynomial, and given the complex value x, this routine improves x by
119  * Laguerre's method until it converges, within the achievable roundoff limit,
120  * to a root of the given polynomial. The number of iterations taken is
121  * returned at its.
122  */
123  void laguer(gmp_complex ** a, int m, gmp_complex * x, int * its, bool type);
124  void computefx(gmp_complex **a, gmp_complex x, int m,
125  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
126  gmp_float &ex, gmp_float &ef);
127  void computegx(gmp_complex **a, gmp_complex x, int m,
128  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
129  gmp_float &ex, gmp_float &ef);
130  void checkimag(gmp_complex *x, gmp_float &e);
131 
132  int var;
133  int tdg;
134 
135  number * coeffs;
136  number * ievpoint;
138 
140 
141  int anz;
143 };
144 //<-
145 
146 class slists; typedef slists * lists;
147 
148 //-> class rootArranger
150 {
151 public:
152  friend lists listOfRoots( rootArranger*, const unsigned int oprec );
153 
154  rootArranger( rootContainer ** _roots,
155  rootContainer ** _mu,
156  const int _howclean = PM_CORRUPT );
158 
159  void solve_all();
160  void arrange();
161 
162  bool success() { return found_roots; }
163 
164 private:
165  rootArranger( const rootArranger & );
166 
169 
170  int howclean;
171  int rc,mc;
173 };
174 
175 
176 
177 //<-
178 
179 //-> simplex computation
180 // (used by sparse matrix construction)
181 #define SIMPLEX_EPS 1.0e-12
182 
183 /** Linear Programming / Linear Optimization using Simplex - Algorithm
184  *
185  * On output, the tableau LiPM is indexed by two arrays of integers.
186  * ipsov[j] contains, for j=1..m, the number i whose original variable
187  * is now represented by row j+1 of LiPM. (left-handed vars in solution)
188  * (first row is the one with the objective function)
189  * izrov[j] contains, for j=1..n, the number i whose original variable
190  * x_i is now a right-handed variable, rep. by column j+1 of LiPM.
191  * These vars are all zero in the solution. The meaning of n<i<n+m1+m2
192  * is the same as above.
193  */
194 class simplex
195 {
196 public:
197 
198  int m; // number of constraints, make sure m == m1 + m2 + m3 !!
199  int n; // # of independent variables
200  int m1,m2,m3; // constraints <=, >= and ==
201  int icase; // == 0: finite solution found;
202  // == +1 objective funtion unbound; == -1: no solution
203  int *izrov,*iposv;
204 
205  mprfloat **LiPM; // the matrix (of size [m+2, n+1])
206 
207  /** #rows should be >= m+2, #cols >= n+1
208  */
209  simplex( int rows, int cols );
210  ~simplex();
211 
214  intvec * posvToIV();
215  intvec * zrovToIV();
216 
217  void compute();
218 
219 private:
220  simplex( const simplex & );
221  void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax );
222  void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 );
223  void simp3( mprfloat **a, int i1, int k1, int ip, int kp );
224 
226 };
227 
228 //<-
229 
230 #endif /*MPR_NUMERIC_H*/
231 
232 // local Variables: ***
233 // folded-file: t ***
234 // compile-command-1: "make installg" ***
235 // compile-command-2: "make install" ***
236 // End: ***
PM_CORRUPT
#define PM_CORRUPT
Definition: mpr_numeric.h:21
rootContainer::coeffs
number * coeffs
Definition: mpr_numeric.h:135
vandermonde::init
void init()
Definition: mpr_numeric.cc:58
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:151
rootContainer::onepoly
@ onepoly
Definition: mpr_numeric.h:68
mpr_complex.h
ip_smatrix
Definition: matpol.h:14
simplex::m
int m
Definition: mpr_numeric.h:198
j
int j
Definition: facHensel.cc:105
simplex::simp3
void simp3(mprfloat **a, int i1, int k1, int ip, int kp)
Definition: mpr_numeric.cc:1341
rootContainer::rootType
rootType
Definition: mpr_numeric.h:68
rootContainer::operator[]
gmp_complex & operator[](const int i)
Definition: mpr_numeric.h:82
rootContainer::getPoly
poly getPoly()
Definition: mpr_numeric.cc:338
k
int k
Definition: cfEzgcd.cc:92
rootContainer::divlin
void divlin(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:642
vandermonde::maxdeg
long maxdeg
Definition: mpr_numeric.h:51
x
Variable x
Definition: cfModGcd.cc:4023
rootContainer::isfloat
bool isfloat(gmp_complex **a)
Definition: mpr_numeric.cc:629
vandermonde::vandermonde
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)
Definition: mpr_numeric.cc:40
simplex::izrov
int * izrov
Definition: mpr_numeric.h:203
rootContainer::solver
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:441
rootContainer::ievpoint
number * ievpoint
Definition: mpr_numeric.h:136
rootContainer::evPointCoord
gmp_complex & evPointCoord(const int i)
Definition: mpr_numeric.cc:392
simplex::LiPM
mprfloat ** LiPM
Definition: mpr_numeric.h:205
vandermonde::x
number * x
Definition: mpr_numeric.h:55
simplex::simp2
void simp2(mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
Definition: mpr_numeric.cc:1302
rootContainer::~rootContainer
~rootContainer()
Definition: mpr_numeric.cc:282
vandermonde
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
rootArranger::mc
int mc
Definition: mpr_numeric.h:171
rootContainer::found_roots
bool found_roots
Definition: mpr_numeric.h:142
rootContainer::rt
rootType rt
Definition: mpr_numeric.h:137
simplex::~simplex
~simplex()
Definition: mpr_numeric.cc:1001
vandermonde::numvec2poly
poly numvec2poly(const number *q)
Definition: mpr_numeric.cc:98
simplex::posvToIV
intvec * posvToIV()
Definition: mpr_numeric.cc:1077
simplex::m1
int m1
Definition: mpr_numeric.h:200
rootContainer::swapRoots
bool swapRoots(const int from, const int to)
Definition: mpr_numeric.cc:421
rootContainer::rootContainer
rootContainer()
Definition: mpr_numeric.cc:269
simplex::simp1
void simp1(mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
Definition: mpr_numeric.cc:1267
i
int i
Definition: cfEzgcd.cc:125
mpr_global.h
rootContainer::getLDim
int getLDim()
Definition: mpr_numeric.h:96
BOOLEAN
int BOOLEAN
Definition: auxiliary.h:85
simplex::m2
int m2
Definition: mpr_numeric.h:200
rootContainer::anz
int anz
Definition: mpr_numeric.h:141
simplex::icase
int icase
Definition: mpr_numeric.h:201
rootContainer::getRoot
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
rootArranger::found_roots
bool found_roots
Definition: mpr_numeric.h:172
simplex::simplex
simplex(int rows, int cols)
#rows should be >= m+2, #cols >= n+1
Definition: mpr_numeric.cc:976
rootContainer::theroots
gmp_complex ** theroots
Definition: mpr_numeric.h:139
rootContainer::tdg
int tdg
Definition: mpr_numeric.h:133
intvec
Definition: intvec.h:17
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:554
rootArranger::rootArranger
rootArranger(rootContainer **_roots, rootContainer **_mu, const int _howclean=PM_CORRUPT)
Definition: mpr_numeric.cc:852
rootContainer::checkimag
void checkimag(gmp_complex *x, gmp_float &e)
Definition: mpr_numeric.cc:621
rootContainer::sortre
void sortre(gmp_complex **r, int l, int u, int inc)
Definition: mpr_numeric.cc:758
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 ad0..tdg of the polynomial this routine succe...
Definition: mpr_numeric.cc:471
rootArranger
Definition: mpr_numeric.h:149
vandermonde::l
long l
Definition: mpr_numeric.h:52
rootArranger::rc
int rc
Definition: mpr_numeric.h:171
rootContainer::fillContainer
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:304
rootArranger::~rootArranger
~rootArranger()
Definition: mpr_numeric.h:157
vandermonde::~vandermonde
~vandermonde()
Definition: mpr_numeric.cc:51
simplex
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
slists
Definition: lists.h:22
simplex::LiPM_rows
int LiPM_rows
Definition: mpr_numeric.h:225
PM_NONE
#define PM_NONE
Definition: mpr_numeric.h:19
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:826
rootArranger::roots
rootContainer ** roots
Definition: mpr_numeric.h:167
rootContainer::getAnzRoots
int getAnzRoots()
Definition: mpr_numeric.h:97
rootContainer::solvequad
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
Definition: mpr_numeric.cc:686
simplex::m3
int m3
Definition: mpr_numeric.h:200
mprfloat
double mprfloat
Definition: mpr_global.h:17
simplex::zrovToIV
intvec * zrovToIV()
Definition: mpr_numeric.cc:1088
rootContainer::cspecialmu
@ cspecialmu
Definition: mpr_numeric.h:68
rootContainer::none
@ none
Definition: mpr_numeric.h:68
rootContainer::divquad
void divquad(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:662
rootArranger::success
bool success()
Definition: mpr_numeric.h:162
rootContainer::det
@ det
Definition: mpr_numeric.h:68
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:805
rootArranger::arrange
void arrange()
Definition: mpr_numeric.cc:887
m
int m
Definition: cfEzgcd.cc:121
lists
slists * lists
Definition: mpr_numeric.h:146
simplex::LiPM_cols
int LiPM_cols
Definition: mpr_numeric.h:225
l
int l
Definition: cfEzgcd.cc:93
simplex::compute
void compute()
Definition: mpr_numeric.cc:1099
rootArranger::solve_all
void solve_all()
Definition: mpr_numeric.cc:862
rootArranger::listOfRoots
friend lists listOfRoots(rootArranger *, const unsigned int oprec)
Definition: ipshell.cc:5003
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
vandermonde::homog
bool homog
Definition: mpr_numeric.h:57
simplex::mapFromMatrix
BOOLEAN mapFromMatrix(matrix m)
Definition: mpr_numeric.cc:1015
rootArranger::mu
rootContainer ** mu
Definition: mpr_numeric.h:168
rootContainer::cspecial
@ cspecial
Definition: mpr_numeric.h:68
rootContainer::var
int var
Definition: mpr_numeric.h:132
rootContainer::getAnzElems
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer::sortroots
void sortroots(gmp_complex **roots, int r, int c, bool isf)
Definition: mpr_numeric.cc:739
rootContainer
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
simplex::n
int n
Definition: mpr_numeric.h:199
rootArranger::howclean
int howclean
Definition: mpr_numeric.h:170
vandermonde::p
number * p
Definition: mpr_numeric.h:54
vandermonde::n
long n
Definition: mpr_numeric.h:49
numbers.h
simplex::mapToMatrix
matrix mapToMatrix(matrix m)
Definition: mpr_numeric.cc:1044
gmp_float
Definition: mpr_complex.h:31
vandermonde::cn
long cn
Definition: mpr_numeric.h:50
simplex::iposv
int * iposv
Definition: mpr_numeric.h:203
gmp_complex
gmp_complex numbers based on
Definition: mpr_complex.h:178