source: git/numeric/mpr_numeric.h @ 6caad65

jengelh-datetimespielwiese
Last change on this file since 6caad65 was 6caad65, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
warnings elimination chg: uninitialized fake variables chg: "'const' type qualifier on return type has no effect" for numeric/ chg: more minor warnings elimination
  • Property mode set to 100644
File size: 6.5 KB
Line 
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 */
28class vandermonde
29{
30public:
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
45private:
46  void init();
47
48private:
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 */
65class rootContainer
66{
67public:
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
99private:
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;
137  rootType rt;
138
139  gmp_complex ** theroots;
140
141  int anz;
142  bool found_roots;
143};
144//<-
145
146class slists; typedef slists * lists;
147
148//-> class rootArranger
149class rootArranger
150{
151public:
152  friend lists listOfRoots( rootArranger*, const unsigned int oprec );
153   
154  rootArranger( rootContainer ** _roots,
155                rootContainer ** _mu,
156                const int _howclean = PM_CORRUPT );
157  ~rootArranger() {}
158
159  void solve_all();
160  void arrange();
161
162  bool success() { return found_roots; }
163
164private:
165  rootArranger( const rootArranger & );
166
167  rootContainer ** roots;
168  rootContainer ** mu;
169
170  int howclean;
171  int rc,mc;
172  bool found_roots;
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 */
194class simplex
195{
196public:
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
219private:
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: ***
Note: See TracBrowser for help on using the repository browser.