source: git/kernel/mpr_numeric.h @ 06662e

spielwiese
Last change on this file since 06662e was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 6.4 KB
Line 
1#ifndef MPR_NUMERIC_H
2#define MPR_NUMERIC_H
3/****************************************
4*  Computer Algebra System SINGULAR     *
5****************************************/
6
7/* $Id$ */
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  {
85    return *theroots[i];
86  }
87  gmp_complex & evPointCoord( const int i );
88
89  inline gmp_complex * getRoot( const int i )
90  {
91    return theroots[i];
92  }
93
94  bool swapRoots( const int from, const int to );
95
96  int getAnzElems() { return anz; }
97  int getLDim() { return anz; }
98  int getAnzRoots() { return tdg; }
99
100private:
101  rootContainer( const rootContainer & v );
102
103  /** Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg]
104   * (generated from the number coefficients coeffs[0..tdg]) of the polynomial
105   * this routine successively calls "laguer" and finds all m complex roots in
106   * roots[0..tdg]. The bool var "polish" should be input as "true" if polishing
107   * (also by "laguer") is desired, "false" if the roots will be subsequently
108   * polished by other means.
109   */
110  bool laguer_driver( gmp_complex ** a, gmp_complex ** roots, bool polish = true );
111  bool isfloat(gmp_complex **a);
112  void divlin(gmp_complex **a, gmp_complex x, int j);
113  void divquad(gmp_complex **a, gmp_complex x, int j);
114  void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j);
115  void sortroots(gmp_complex **roots, int r, int c, bool isf);
116  void sortre(gmp_complex **r, int l, int u, int inc);
117
118  /** Given the degree m and the m+1 complex coefficients a[0..m] of the
119   * polynomial, and given the complex value x, this routine improves x by
120   * Laguerre's method until it converges, within the achievable roundoff limit,
121   * to a root of the given polynomial. The number of iterations taken is
122   * returned at its.
123   */
124  void laguer(gmp_complex ** a, int m, gmp_complex * x, int * its, bool type);
125  void computefx(gmp_complex **a, gmp_complex x, int m,
126                gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
127                gmp_float &ex, gmp_float &ef);
128  void computegx(gmp_complex **a, gmp_complex x, int m,
129                gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
130                gmp_float &ex, gmp_float &ef);
131  void checkimag(gmp_complex *x, gmp_float &e);
132
133  int var;
134  int tdg;
135
136  number * coeffs;
137  number * ievpoint;
138  rootType rt;
139
140  gmp_complex ** theroots;
141
142  int anz;
143  bool found_roots;
144};
145//<-
146
147//-> class rootArranger
148class rootArranger
149{
150public:
151  rootArranger( rootContainer ** _roots,
152                rootContainer ** _mu,
153                const int _howclean = PM_CORRUPT );
154  ~rootArranger() {}
155
156  void solve_all();
157  void arrange();
158
159  lists listOfRoots( const unsigned int oprec );
160
161  const bool success() { return found_roots; }
162
163private:
164  rootArranger( const rootArranger & );
165
166  rootContainer ** roots;
167  rootContainer ** mu;
168
169  int howclean;
170  int rc,mc;
171  bool found_roots;
172};
173//<-
174
175//-> simplex computation
176//   (used by sparse matrix construction)
177#define SIMPLEX_EPS 1.0e-12
178
179/** Linear Programming / Linear Optimization using Simplex - Algorithm
180 *
181 * On output, the tableau LiPM is indexed by two arrays of integers.
182 * ipsov[j] contains, for j=1..m, the number i whose original variable
183 * is now represented by row j+1 of LiPM. (left-handed vars in solution)
184 * (first row is the one with the objective function)
185 * izrov[j] contains, for j=1..n, the number i whose original variable
186 * x_i is now a right-handed variable, rep. by column j+1 of LiPM.
187 * These vars are all zero in the solution. The meaning of n<i<n+m1+m2
188 * is the same as above.
189 */
190class simplex
191{
192public:
193
194  int m;         // number of constraints, make sure m == m1 + m2 + m3 !!
195  int n;         // # of independent variables
196  int m1,m2,m3;  // constraints <=, >= and ==
197  int icase;     // == 0: finite solution found;
198                 // == +1 objective funtion unbound; == -1: no solution
199  int *izrov,*iposv;
200
201  mprfloat **LiPM; // the matrix (of size [m+2, n+1])
202
203  /** #rows should be >= m+2, #cols >= n+1
204   */
205  simplex( int rows, int cols );
206  ~simplex();
207
208  BOOLEAN mapFromMatrix( matrix m );
209  matrix mapToMatrix( matrix m );
210  intvec * posvToIV();
211  intvec * zrovToIV();
212
213  void compute();
214
215private:
216  simplex( const simplex & );
217  void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax );
218  void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 );
219  void simp3( mprfloat **a, int i1, int k1, int ip, int kp );
220
221  int LiPM_cols,LiPM_rows;
222};
223
224//<-
225
226#endif /*MPR_NUMERIC_H*/
227
228// local Variables: ***
229// folded-file: t ***
230// compile-command-1: "make installg" ***
231// compile-command-2: "make install" ***
232// End: ***
Note: See TracBrowser for help on using the repository browser.