My Project
Loading...
Searching...
No Matches
mpr_numeric.h
Go to the documentation of this file.
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 */
29{
30public:
31 vandermonde( const long _cn, const long _n,
32 const long _maxdeg, number *_p, const bool _homog = true );
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 */
66{
67public:
69
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:
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;
138
140
141 int anz;
143};
144//<-
145
146class slists; typedef slists * lists;
147
148//-> class 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 );
158
159 void solve_all();
160 void arrange();
161
162 bool success() { return found_roots; }
163
164private:
166
169
171 int rc,mc;
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 */
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 function unbound; == -1: no solution
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
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
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: ***
int BOOLEAN
Definition: auxiliary.h:87
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
gmp_complex numbers based on
Definition: mpr_complex.h:179
Definition: intvec.h:23
void solve_all()
Definition: mpr_numeric.cc:858
friend lists listOfRoots(rootArranger *, const unsigned int oprec)
Definition: ipshell.cc:5079
rootContainer ** roots
Definition: mpr_numeric.h:167
rootArranger(const rootArranger &)
rootContainer ** mu
Definition: mpr_numeric.h:168
bool found_roots
Definition: mpr_numeric.h:172
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:883
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
void sortre(gmp_complex **r, int l, int u, int inc)
Definition: mpr_numeric.cc:754
gmp_complex & operator[](const int i)
Definition: mpr_numeric.h:82
bool laguer_driver(gmp_complex **a, gmp_complex **roots, bool polish=true)
Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg] (generated from the number coeffic...
Definition: mpr_numeric.cc:467
void computegx(gmp_complex **a, gmp_complex x, int m, gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, gmp_float &ex, gmp_float &ef)
Definition: mpr_numeric.cc:822
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:300
void laguer(gmp_complex **a, int m, gmp_complex *x, int *its, bool type)
Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial, and given the complex ...
Definition: mpr_numeric.cc:550
void divlin(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:638
void sortroots(gmp_complex **roots, int r, int c, bool isf)
Definition: mpr_numeric.cc:735
rootType rt
Definition: mpr_numeric.h:137
void divquad(gmp_complex **a, gmp_complex x, int j)
Definition: mpr_numeric.cc:658
bool swapRoots(const int from, const int to)
Definition: mpr_numeric.cc:417
gmp_complex ** theroots
Definition: mpr_numeric.h:139
void computefx(gmp_complex **a, gmp_complex x, int m, gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, gmp_float &ex, gmp_float &ef)
Definition: mpr_numeric.cc:801
number * coeffs
Definition: mpr_numeric.h:135
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
Definition: mpr_numeric.cc:682
gmp_complex & evPointCoord(const int i)
Definition: mpr_numeric.cc:388
bool isfloat(gmp_complex **a)
Definition: mpr_numeric.cc:625
rootContainer(const rootContainer &v)
void checkimag(gmp_complex *x, gmp_float &e)
Definition: mpr_numeric.cc:617
int getAnzRoots()
Definition: mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:437
number * ievpoint
Definition: mpr_numeric.h:136
int getAnzElems()
Definition: mpr_numeric.h:95
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:195
intvec * zrovToIV()
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int * iposv
Definition: mpr_numeric.h:203
int LiPM_rows
Definition: mpr_numeric.h:225
BOOLEAN mapFromMatrix(matrix m)
int * izrov
Definition: mpr_numeric.h:203
int icase
Definition: mpr_numeric.h:201
simplex(const simplex &)
void compute()
int LiPM_cols
Definition: mpr_numeric.h:225
void simp2(mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
void simp3(mprfloat **a, int i1, int k1, int ip, int kp)
void simp1(mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
matrix mapToMatrix(matrix m)
intvec * posvToIV()
Definition: lists.h:24
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
void init()
Definition: mpr_numeric.cc:53
poly numvec2poly(const number *q)
Definition: mpr_numeric.cc:93
number * x
Definition: mpr_numeric.h:55
number * p
Definition: mpr_numeric.h:54
number * interpolateDense(const number *q)
Solves the Vandermode linear system \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
Definition: mpr_numeric.cc:146
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
double mprfloat
Definition: mpr_global.h:17
#define PM_NONE
Definition: mpr_numeric.h:19
#define PM_CORRUPT
Definition: mpr_numeric.h:21
slists * lists
Definition: mpr_numeric.h:146