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.1.1.1 2003-10-06 12:15:58 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 | 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 | |
---|
98 | private: |
---|
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 |
---|
146 | class rootArranger |
---|
147 | { |
---|
148 | public: |
---|
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 | |
---|
161 | private: |
---|
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 | */ |
---|
188 | class simplex |
---|
189 | { |
---|
190 | public: |
---|
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 | |
---|
213 | private: |
---|
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: *** |
---|