source:git/kernel/mpr_numeric.h@91f1a3

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