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.3 2007-05-24 17:46:04 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 | */ |
---|
29 | class vandermonde |
---|
30 | { |
---|
31 | public: |
---|
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 | |
---|
46 | private: |
---|
47 | void init(); |
---|
48 | |
---|
49 | private: |
---|
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 | */ |
---|
66 | class rootContainer |
---|
67 | { |
---|
68 | public: |
---|
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 | |
---|
100 | private: |
---|
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 |
---|
148 | class rootArranger |
---|
149 | { |
---|
150 | public: |
---|
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 | |
---|
163 | private: |
---|
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 | */ |
---|
190 | class simplex |
---|
191 | { |
---|
192 | public: |
---|
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 | |
---|
215 | private: |
---|
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: *** |
---|