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:
68  enum rootType { none, cspecial, cspecialmu, det, onepoly };
69 
70  rootContainer();
71  ~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 
212  BOOLEAN mapFromMatrix( matrix m );
213  matrix mapToMatrix( matrix m );
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 
225  int LiPM_cols,LiPM_rows;
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: ***
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
int j
Definition: facHensel.cc:105
number * interpolateDense(const number *q)
Solves the Vandermode linear system {i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
Definition: mpr_numeric.cc:151
rootContainer ** mu
Definition: mpr_numeric.h:168
#define PM_CORRUPT
Definition: mpr_numeric.h:21
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
number * ievpoint
Definition: mpr_numeric.h:136
Definition: lists.h:22
number * coeffs
Definition: mpr_numeric.h:135
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
void init()
Definition: mpr_numeric.cc:58
rootType rt
Definition: mpr_numeric.h:137
int k
Definition: cfEzgcd.cc:92
gmp_complex numbers based on
Definition: mpr_complex.h:178
int getAnzElems()
Definition: mpr_numeric.h:95
poly numvec2poly(const number *q)
Definition: mpr_numeric.cc:98
double mprfloat
Definition: mpr_global.h:17
number * x
Definition: mpr_numeric.h:55
Definition: intvec.h:17
#define PM_NONE
Definition: mpr_numeric.h:19
bool success()
Definition: mpr_numeric.h:162
int m
Definition: cfEzgcd.cc:121
int i
Definition: cfEzgcd.cc:125
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
bool found_roots
Definition: mpr_numeric.h:172
gmp_complex ** theroots
Definition: mpr_numeric.h:139
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
mprfloat ** LiPM
Definition: mpr_numeric.h:205
slists * lists
Definition: mpr_numeric.h:146
int getAnzRoots()
Definition: mpr_numeric.h:97
int LiPM_rows
Definition: mpr_numeric.h:225
number * p
Definition: mpr_numeric.h:54
rootContainer ** roots
Definition: mpr_numeric.h:167
int icase
Definition: mpr_numeric.h:201
int * izrov
Definition: mpr_numeric.h:203
int BOOLEAN
Definition: auxiliary.h:85
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:5003
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)
Definition: mpr_numeric.cc:40