1 | /**************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | ****************************************/ |
---|
4 | |
---|
5 | /* $Id: mpr_numeric.cc,v 1.3 2005-05-04 07:23:17 Singular Exp $ */ |
---|
6 | |
---|
7 | /* |
---|
8 | * ABSTRACT - multipolynomial resultants - numeric stuff |
---|
9 | * ( root finder, vandermonde system solver, simplex ) |
---|
10 | */ |
---|
11 | |
---|
12 | #include "mod2.h" |
---|
13 | //#ifdef HAVE_MPR |
---|
14 | |
---|
15 | //#define mprDEBUG_ALL |
---|
16 | |
---|
17 | //-> includes |
---|
18 | #include "structs.h" |
---|
19 | #include "febase.h" |
---|
20 | #include "omalloc.h" |
---|
21 | #include "numbers.h" |
---|
22 | #include "polys.h" |
---|
23 | #include "ideals.h" |
---|
24 | #include "intvec.h" |
---|
25 | #include "matpol.h" |
---|
26 | #include "ring.h" |
---|
27 | //#include "longrat.h" |
---|
28 | |
---|
29 | #include <math.h> |
---|
30 | #include "mpr_numeric.h" |
---|
31 | |
---|
32 | extern size_t gmp_output_digits; |
---|
33 | //<- |
---|
34 | |
---|
35 | extern void nPrint(number n); // for debugging output |
---|
36 | |
---|
37 | //----------------------------------------------------------------------------- |
---|
38 | //-------------- vandermonde system solver ------------------------------------ |
---|
39 | //----------------------------------------------------------------------------- |
---|
40 | |
---|
41 | //-> vandermonde::* |
---|
42 | vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg, |
---|
43 | number *_p, const bool _homog ) |
---|
44 | : n(_n), cn(_cn), maxdeg(_maxdeg), p(_p), homog(_homog) |
---|
45 | { |
---|
46 | long j; |
---|
47 | l= (long)pow((double)maxdeg+1,(int)n); |
---|
48 | x= (number *)omAlloc( cn * sizeof(number) ); |
---|
49 | for ( j= 0; j < cn; j++ ) x[j]= nInit(1); |
---|
50 | init(); |
---|
51 | } |
---|
52 | |
---|
53 | vandermonde::~vandermonde() |
---|
54 | { |
---|
55 | int j; |
---|
56 | for ( j= 0; j < cn; j++ ) nDelete( x+j ); |
---|
57 | omFreeSize( (ADDRESS)x, cn * sizeof( number ) ); |
---|
58 | } |
---|
59 | |
---|
60 | void vandermonde::init() |
---|
61 | { |
---|
62 | int j; |
---|
63 | long i,c,sum; |
---|
64 | number tmp,tmp1; |
---|
65 | |
---|
66 | c=0; |
---|
67 | sum=0; |
---|
68 | |
---|
69 | intvec exp( n ); |
---|
70 | for ( j= 0; j < n; j++ ) exp[j]=0; |
---|
71 | |
---|
72 | for ( i= 0; i < l; i++ ) |
---|
73 | { |
---|
74 | if ( !homog || (sum == maxdeg) ) |
---|
75 | { |
---|
76 | for ( j= 0; j < n; j++ ) |
---|
77 | { |
---|
78 | nPower( p[j], exp[j], &tmp ); |
---|
79 | tmp1 = nMult( tmp, x[c] ); |
---|
80 | x[c]= tmp1; |
---|
81 | nDelete( &tmp ); |
---|
82 | } |
---|
83 | c++; |
---|
84 | } |
---|
85 | exp[0]++; |
---|
86 | sum=0; |
---|
87 | for ( j= 0; j < n - 1; j++ ) |
---|
88 | { |
---|
89 | if ( exp[j] > maxdeg ) |
---|
90 | { |
---|
91 | exp[j]= 0; |
---|
92 | exp[j + 1]++; |
---|
93 | } |
---|
94 | sum+= exp[j]; |
---|
95 | } |
---|
96 | sum+= exp[n - 1]; |
---|
97 | } |
---|
98 | } |
---|
99 | |
---|
100 | poly vandermonde::numvec2poly( const number * q ) |
---|
101 | { |
---|
102 | int j; |
---|
103 | long i,c,sum; |
---|
104 | number tmp; |
---|
105 | |
---|
106 | poly pnew,pit=NULL; |
---|
107 | |
---|
108 | c=0; |
---|
109 | sum=0; |
---|
110 | |
---|
111 | int *exp= (int *) omAlloc( (n+1) * sizeof(int) ); |
---|
112 | |
---|
113 | for ( j= 0; j < n+1; j++ ) exp[j]=0; |
---|
114 | |
---|
115 | for ( i= 0; i < l; i++ ) |
---|
116 | { |
---|
117 | if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) ) |
---|
118 | { |
---|
119 | pnew= pOne(); |
---|
120 | pSetCoeff(pnew,q[i]); |
---|
121 | pSetExpV(pnew,exp); |
---|
122 | if ( pit ) |
---|
123 | { |
---|
124 | pNext(pnew)= pit; |
---|
125 | pit= pnew; |
---|
126 | } |
---|
127 | else |
---|
128 | { |
---|
129 | pit= pnew; |
---|
130 | pNext(pnew)= NULL; |
---|
131 | } |
---|
132 | pSetm(pit); |
---|
133 | } |
---|
134 | exp[1]++; |
---|
135 | sum=0; |
---|
136 | for ( j= 1; j < n; j++ ) |
---|
137 | { |
---|
138 | if ( exp[j] > maxdeg ) |
---|
139 | { |
---|
140 | exp[j]= 0; |
---|
141 | exp[j + 1]++; |
---|
142 | } |
---|
143 | sum+= exp[j]; |
---|
144 | } |
---|
145 | sum+= exp[n]; |
---|
146 | } |
---|
147 | |
---|
148 | omFreeSize( (ADDRESS) exp, (n+1) * sizeof(int) ); |
---|
149 | |
---|
150 | pSortAdd(pit); |
---|
151 | return pit; |
---|
152 | } |
---|
153 | |
---|
154 | number * vandermonde::interpolateDense( const number * q ) |
---|
155 | { |
---|
156 | int i,j,k; |
---|
157 | number newnum,tmp1; |
---|
158 | number b,t,xx,s; |
---|
159 | number *c; |
---|
160 | number *w; |
---|
161 | |
---|
162 | b=t=xx=s=tmp1=NULL; |
---|
163 | |
---|
164 | w= (number *)omAlloc( cn * sizeof(number) ); |
---|
165 | c= (number *)omAlloc( cn * sizeof(number) ); |
---|
166 | for ( j= 0; j < cn; j++ ) |
---|
167 | { |
---|
168 | w[j]= nInit(0); |
---|
169 | c[j]= nInit(0); |
---|
170 | } |
---|
171 | |
---|
172 | if ( cn == 1 ) |
---|
173 | { |
---|
174 | nDelete( &w[0] ); |
---|
175 | w[0]= nCopy(q[0]); |
---|
176 | } |
---|
177 | else |
---|
178 | { |
---|
179 | nDelete( &c[cn-1] ); |
---|
180 | c[cn-1]= nCopy(x[0]); |
---|
181 | c[cn-1]= nNeg(c[cn-1]); // c[cn]= -x[1] |
---|
182 | |
---|
183 | for ( i= 1; i < cn; i++ ) { // i=2; i <= cn |
---|
184 | nDelete( &xx ); |
---|
185 | xx= nCopy(x[i]); |
---|
186 | xx= nNeg(xx); // xx= -x[i] |
---|
187 | |
---|
188 | for ( j= (cn-i-1); j <= (cn-2); j++) { // j=(cn+1-i); j <= (cn-1) |
---|
189 | nDelete( &tmp1 ); |
---|
190 | tmp1= nMult( xx, c[j+1] ); // c[j]= c[j] + (xx * c[j+1]) |
---|
191 | newnum= nAdd( c[j], tmp1 ); |
---|
192 | nDelete( &c[j] ); |
---|
193 | c[j]= newnum; |
---|
194 | } |
---|
195 | |
---|
196 | newnum= nAdd( xx, c[cn-1] ); // c[cn-1]= c[cn-1] + xx |
---|
197 | nDelete( &c[cn-1] ); |
---|
198 | c[cn-1]= newnum; |
---|
199 | } |
---|
200 | |
---|
201 | for ( i= 0; i < cn; i++ ) { // i=1; i <= cn |
---|
202 | nDelete( &xx ); |
---|
203 | xx= nCopy(x[i]); // xx= x[i] |
---|
204 | |
---|
205 | nDelete( &t ); |
---|
206 | t= nInit( 1 ); // t= b= 1 |
---|
207 | nDelete( &b ); |
---|
208 | b= nInit( 1 ); |
---|
209 | nDelete( &s ); // s= q[cn-1] |
---|
210 | s= nCopy( q[cn-1] ); |
---|
211 | |
---|
212 | for ( k= cn-1; k >= 1; k-- ) { // k=cn; k >= 2 |
---|
213 | nDelete( &tmp1 ); |
---|
214 | tmp1= nMult( xx, b ); // b= c[k] + (xx * b) |
---|
215 | nDelete( &b ); |
---|
216 | b= nAdd( c[k], tmp1 ); |
---|
217 | |
---|
218 | nDelete( &tmp1 ); |
---|
219 | tmp1= nMult( q[k-1], b ); // s= s + (q[k-1] * b) |
---|
220 | newnum= nAdd( s, tmp1 ); |
---|
221 | nDelete( &s ); |
---|
222 | s= newnum; |
---|
223 | |
---|
224 | nDelete( &tmp1 ); |
---|
225 | tmp1= nMult( xx, t ); // t= (t * xx) + b |
---|
226 | newnum= nAdd( tmp1, b ); |
---|
227 | nDelete( &t ); |
---|
228 | t= newnum; |
---|
229 | } |
---|
230 | |
---|
231 | if (!nIsZero(t)) |
---|
232 | { |
---|
233 | nDelete( &w[i] ); // w[i]= s/t |
---|
234 | w[i]= nDiv( s, t ); |
---|
235 | nNormalize( w[i] ); |
---|
236 | } |
---|
237 | |
---|
238 | mprSTICKYPROT(ST_VANDER_STEP); |
---|
239 | } |
---|
240 | } |
---|
241 | mprSTICKYPROT("\n"); |
---|
242 | |
---|
243 | // free mem |
---|
244 | for ( j= 0; j < cn; j++ ) nDelete( c+j ); |
---|
245 | omFreeSize( (ADDRESS)c, cn * sizeof( number ) ); |
---|
246 | |
---|
247 | nDelete( &tmp1 ); |
---|
248 | nDelete( &s ); |
---|
249 | nDelete( &t ); |
---|
250 | nDelete( &b ); |
---|
251 | nDelete( &xx ); |
---|
252 | |
---|
253 | // makes quotiens smaller |
---|
254 | for ( j= 0; j < cn; j++ ) nNormalize( w[j] ); |
---|
255 | |
---|
256 | return w; |
---|
257 | } |
---|
258 | //<- |
---|
259 | |
---|
260 | //----------------------------------------------------------------------------- |
---|
261 | //-------------- rootContainer ------------------------------------------------ |
---|
262 | //----------------------------------------------------------------------------- |
---|
263 | |
---|
264 | //-> definitions |
---|
265 | #define MR 8 // never change this value |
---|
266 | #define MT 5 |
---|
267 | #define MMOD (MT*MR) |
---|
268 | #define MAXIT (5*MMOD) // max number of iterations in laguer root finder |
---|
269 | |
---|
270 | |
---|
271 | //-> rootContainer::rootContainer() |
---|
272 | rootContainer::rootContainer() |
---|
273 | { |
---|
274 | rt=none; |
---|
275 | |
---|
276 | coeffs= NULL; |
---|
277 | ievpoint= NULL; |
---|
278 | theroots= NULL; |
---|
279 | |
---|
280 | found_roots= false; |
---|
281 | } |
---|
282 | //<- |
---|
283 | |
---|
284 | //-> rootContainer::~rootContainer() |
---|
285 | rootContainer::~rootContainer() |
---|
286 | { |
---|
287 | int i; |
---|
288 | int n= pVariables; |
---|
289 | |
---|
290 | // free coeffs, ievpoint |
---|
291 | if ( ievpoint != NULL ) |
---|
292 | { |
---|
293 | for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i ); |
---|
294 | omFreeSize( (ADDRESS)ievpoint, (anz+2) * sizeof( number ) ); |
---|
295 | } |
---|
296 | |
---|
297 | for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i ); |
---|
298 | omFreeSize( (ADDRESS)coeffs, (tdg+1) * sizeof( number ) ); |
---|
299 | |
---|
300 | // theroots löschen |
---|
301 | for ( i=0; i < tdg; i++ ) delete theroots[i]; |
---|
302 | omFreeSize( (ADDRESS) theroots, (tdg)*sizeof(gmp_complex*) ); |
---|
303 | |
---|
304 | //mprPROTnl("~rootContainer()"); |
---|
305 | } |
---|
306 | //<- |
---|
307 | |
---|
308 | //-> void rootContainer::fillContainer( ... ) |
---|
309 | void rootContainer::fillContainer( number *_coeffs, number *_ievpoint, |
---|
310 | const int _var, const int _tdg, |
---|
311 | const rootType _rt, const int _anz ) |
---|
312 | { |
---|
313 | int i; |
---|
314 | number nn= nInit(0); |
---|
315 | var=_var; |
---|
316 | tdg=_tdg; |
---|
317 | coeffs=_coeffs; |
---|
318 | rt=_rt; |
---|
319 | anz=_anz; |
---|
320 | |
---|
321 | for ( i=0; i <= tdg; i++ ) |
---|
322 | { |
---|
323 | if ( nEqual(coeffs[i],nn) ) |
---|
324 | { |
---|
325 | nDelete( &coeffs[i] ); |
---|
326 | coeffs[i]=NULL; |
---|
327 | } |
---|
328 | } |
---|
329 | nDelete( &nn ); |
---|
330 | |
---|
331 | if ( rt == cspecialmu && _ievpoint ) { // copy ievpoint |
---|
332 | ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) ); |
---|
333 | for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] ); |
---|
334 | } |
---|
335 | |
---|
336 | theroots= NULL; |
---|
337 | found_roots= false; |
---|
338 | } |
---|
339 | //<- |
---|
340 | |
---|
341 | //-> poly rootContainer::getPoly() |
---|
342 | poly rootContainer::getPoly() |
---|
343 | { |
---|
344 | int i; |
---|
345 | |
---|
346 | poly result= NULL; |
---|
347 | poly ppos; |
---|
348 | |
---|
349 | if ( (rt == cspecial) || ( rt == cspecialmu ) ) |
---|
350 | { |
---|
351 | for ( i= tdg; i >= 0; i-- ) |
---|
352 | { |
---|
353 | if ( coeffs[i] ) |
---|
354 | { |
---|
355 | poly p= pOne(); |
---|
356 | //pSetExp( p, var+1, i); |
---|
357 | pSetExp( p, 1, i); |
---|
358 | pSetCoeff( p, nCopy( coeffs[i] ) ); |
---|
359 | pSetm( p ); |
---|
360 | if (result) |
---|
361 | { |
---|
362 | ppos->next=p; |
---|
363 | ppos=ppos->next; |
---|
364 | } |
---|
365 | else |
---|
366 | { |
---|
367 | result=p; |
---|
368 | ppos=p; |
---|
369 | } |
---|
370 | |
---|
371 | } |
---|
372 | } |
---|
373 | if (result!=NULL) pSetm( result ); |
---|
374 | } |
---|
375 | |
---|
376 | return result; |
---|
377 | } |
---|
378 | //<- |
---|
379 | |
---|
380 | //-> const gmp_complex & rootContainer::opterator[] ( const int i ) |
---|
381 | // this is now inline! |
---|
382 | // gmp_complex & rootContainer::operator[] ( const int i ) |
---|
383 | // { |
---|
384 | // if ( found_roots && ( i >= 0) && ( i < tdg ) ) |
---|
385 | // { |
---|
386 | // return *theroots[i]; |
---|
387 | // } |
---|
388 | // // warning |
---|
389 | // Warn("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false"); |
---|
390 | // gmp_complex *tmp= new gmp_complex(); |
---|
391 | // return *tmp; |
---|
392 | // } |
---|
393 | //<- |
---|
394 | |
---|
395 | //-> gmp_complex & rootContainer::evPointCoord( int i ) |
---|
396 | gmp_complex & rootContainer::evPointCoord( const int i ) |
---|
397 | { |
---|
398 | if (! ((i >= 0) && (i < anz+2) ) ) |
---|
399 | WarnS("rootContainer::evPointCoord: index out of range"); |
---|
400 | if (ievpoint == NULL) |
---|
401 | WarnS("rootContainer::evPointCoord: ievpoint == NULL"); |
---|
402 | |
---|
403 | if ( (rt == cspecialmu) && found_roots ) { // FIX ME |
---|
404 | if ( ievpoint[i] != NULL ) |
---|
405 | { |
---|
406 | gmp_complex *tmp= new gmp_complex(); |
---|
407 | *tmp= numberToComplex(ievpoint[i]); |
---|
408 | return *tmp; |
---|
409 | } |
---|
410 | else |
---|
411 | { |
---|
412 | Warn("rootContainer::evPointCoord: NULL index %d",i); |
---|
413 | } |
---|
414 | } |
---|
415 | |
---|
416 | // warning |
---|
417 | Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false"); |
---|
418 | gmp_complex *tmp= new gmp_complex(); |
---|
419 | return *tmp; |
---|
420 | } |
---|
421 | //<- |
---|
422 | |
---|
423 | //-> bool rootContainer::changeRoots( int from, int to ) |
---|
424 | bool rootContainer::swapRoots( const int from, const int to ) |
---|
425 | { |
---|
426 | if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) ) |
---|
427 | { |
---|
428 | if ( to != from ) |
---|
429 | { |
---|
430 | gmp_complex tmp( *theroots[from] ); |
---|
431 | *theroots[from]= *theroots[to]; |
---|
432 | *theroots[to]= tmp; |
---|
433 | } |
---|
434 | return true; |
---|
435 | } |
---|
436 | |
---|
437 | // warning |
---|
438 | Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to); |
---|
439 | return false; |
---|
440 | } |
---|
441 | //<- |
---|
442 | |
---|
443 | //-> void rootContainer::solver() |
---|
444 | bool rootContainer::solver( const int polishmode ) |
---|
445 | { |
---|
446 | int i; |
---|
447 | |
---|
448 | // there are maximal tdg roots, so *roots ranges form 0 to tdg-1. |
---|
449 | theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) ); |
---|
450 | for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex(); |
---|
451 | |
---|
452 | // copy the coefficients of type number to type gmp_complex |
---|
453 | gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) ); |
---|
454 | for ( i=0; i <= tdg; i++ ) |
---|
455 | { |
---|
456 | ad[i]= new gmp_complex(); |
---|
457 | if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i] ); |
---|
458 | } |
---|
459 | |
---|
460 | // now solve |
---|
461 | found_roots= laguer_driver( ad, theroots, polishmode != 0 ); |
---|
462 | if (!found_roots) |
---|
463 | WarnS("rootContainer::solver: No roots found!"); |
---|
464 | |
---|
465 | // free memory |
---|
466 | for ( i=0; i <= tdg; i++ ) delete ad[i]; |
---|
467 | omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) ); |
---|
468 | |
---|
469 | return found_roots; |
---|
470 | } |
---|
471 | //<- |
---|
472 | |
---|
473 | //-> gmp_complex* rootContainer::laguer_driver( bool polish ) |
---|
474 | bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish ) |
---|
475 | { |
---|
476 | int i,j,k,its; |
---|
477 | gmp_float zero(0.0); |
---|
478 | gmp_complex x(0.0),o(1.0); |
---|
479 | bool ret= true, isf=isfloat(a), type=true; |
---|
480 | |
---|
481 | gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) ); |
---|
482 | for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] ); |
---|
483 | |
---|
484 | k = 0; |
---|
485 | i = tdg; |
---|
486 | j = i-1; |
---|
487 | while (i>2) |
---|
488 | { |
---|
489 | // run laguer alg |
---|
490 | x = zero; |
---|
491 | laguer(ad, i, &x, &its, type); |
---|
492 | if ( its > MAXIT ) |
---|
493 | { |
---|
494 | type = !type; |
---|
495 | x = zero; |
---|
496 | laguer(ad, i, &x, &its, type); |
---|
497 | } |
---|
498 | |
---|
499 | mprSTICKYPROT(ST_ROOTS_LGSTEP); |
---|
500 | if ( its > MAXIT ) |
---|
501 | { // error |
---|
502 | WarnS("rootContainer::laguer_driver: To many iterations!"); |
---|
503 | ret= false; |
---|
504 | goto theend; |
---|
505 | } |
---|
506 | if ( polish ) |
---|
507 | { |
---|
508 | laguer( a, tdg, &x, &its , type); |
---|
509 | if ( its > MAXIT ) |
---|
510 | { // error |
---|
511 | WarnS("rootContainer::laguer_driver: To many iterations in polish!"); |
---|
512 | ret= false; |
---|
513 | goto theend; |
---|
514 | } |
---|
515 | } |
---|
516 | if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x; |
---|
517 | if (x.imag() == zero) |
---|
518 | { |
---|
519 | *roots[k] = x; |
---|
520 | k++; |
---|
521 | divlin(ad,x,i); |
---|
522 | i--; |
---|
523 | } |
---|
524 | else |
---|
525 | { |
---|
526 | if(isf) |
---|
527 | { |
---|
528 | *roots[j] = x; |
---|
529 | *roots[j-1]= gmp_complex(x.real(),-x.imag()); |
---|
530 | j -= 2; |
---|
531 | divquad(ad,x,i); |
---|
532 | i -= 2; |
---|
533 | } |
---|
534 | else |
---|
535 | { |
---|
536 | *roots[j] = x; |
---|
537 | j--; |
---|
538 | divlin(ad,x,i); |
---|
539 | i--; |
---|
540 | } |
---|
541 | } |
---|
542 | type = !type; |
---|
543 | } |
---|
544 | solvequad(ad,roots,k,j); |
---|
545 | sortroots(roots,k,j,isf); |
---|
546 | |
---|
547 | theend: |
---|
548 | mprSTICKYPROT("\n"); |
---|
549 | for ( i=0; i <= tdg; i++ ) delete ad[i]; |
---|
550 | omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* )); |
---|
551 | |
---|
552 | return ret; |
---|
553 | } |
---|
554 | //<- |
---|
555 | |
---|
556 | //-> void rootContainer::laguer(...) |
---|
557 | void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its, bool type) |
---|
558 | { |
---|
559 | int iter,j; |
---|
560 | gmp_float zero(0.0),one(1.0),deg(m); |
---|
561 | gmp_float abx_g, err_g, fabs; |
---|
562 | gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2; |
---|
563 | gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 }; |
---|
564 | |
---|
565 | gmp_float epss(0.1); |
---|
566 | mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),getGMPFloatDigits()); |
---|
567 | |
---|
568 | for ( iter= 1; iter <= MAXIT; iter++ ) |
---|
569 | { |
---|
570 | mprSTICKYPROT(ST_ROOTS_LG); |
---|
571 | *its=iter; |
---|
572 | if (type) |
---|
573 | computefx(a,*x,m,b,d,f,abx_g,err_g); |
---|
574 | else |
---|
575 | computegx(a,*x,m,b,d,f,abx_g,err_g); |
---|
576 | err_g *= epss; // EPSS; |
---|
577 | |
---|
578 | fabs = abs(b); |
---|
579 | if (fabs <= err_g) |
---|
580 | { |
---|
581 | if ((fabs==zero) || (abs(d)==zero)) return; |
---|
582 | *x -= (b/d); // a last newton-step |
---|
583 | goto ende; |
---|
584 | } |
---|
585 | |
---|
586 | g= d / b; |
---|
587 | g2 = g * g; |
---|
588 | h= g2 - (((f+f) / b )); |
---|
589 | sq= sqrt(( ( h * deg ) - g2 ) * (deg - one)); |
---|
590 | gp= g + sq; |
---|
591 | gm= g - sq; |
---|
592 | if (abs(gp)<abs(gm)) |
---|
593 | { |
---|
594 | dx = deg/gm; |
---|
595 | } |
---|
596 | else |
---|
597 | { |
---|
598 | if((gp.real()==zero)&&(gp.imag()==zero)) |
---|
599 | { |
---|
600 | dx.real(cos((mprfloat)iter)); |
---|
601 | dx.imag(sin((mprfloat)iter)); |
---|
602 | dx = dx*(one+abx_g); |
---|
603 | } |
---|
604 | else |
---|
605 | { |
---|
606 | dx = deg/gp; |
---|
607 | } |
---|
608 | } |
---|
609 | x1= *x - dx; |
---|
610 | |
---|
611 | if (*x == x1) goto ende; |
---|
612 | |
---|
613 | j = iter%MMOD; |
---|
614 | if (j==0) j=MT; |
---|
615 | if ( j % MT ) *x= x1; |
---|
616 | else *x -= ( dx * frac_g[ j / MT ] ); |
---|
617 | } |
---|
618 | |
---|
619 | *its= MAXIT+1; |
---|
620 | ende: |
---|
621 | checkimag(x,epss); |
---|
622 | } |
---|
623 | |
---|
624 | void rootContainer::checkimag(gmp_complex *x, gmp_float &e) |
---|
625 | { |
---|
626 | if(abs(x->imag())<abs(x->real())*e) |
---|
627 | { |
---|
628 | x->imag(0.0); |
---|
629 | } |
---|
630 | } |
---|
631 | |
---|
632 | bool rootContainer::isfloat(gmp_complex **a) |
---|
633 | { |
---|
634 | gmp_float z(0.0); |
---|
635 | gmp_complex *b; |
---|
636 | for (int i=tdg; i >= 0; i-- ) |
---|
637 | { |
---|
638 | b = &(*a[i]); |
---|
639 | if (!(b->imag()==z)) |
---|
640 | return false; |
---|
641 | } |
---|
642 | return true; |
---|
643 | } |
---|
644 | |
---|
645 | void rootContainer::divlin(gmp_complex **a, gmp_complex x, int j) |
---|
646 | { |
---|
647 | int i; |
---|
648 | gmp_float o(1.0); |
---|
649 | |
---|
650 | if (abs(x)<o) |
---|
651 | { |
---|
652 | for (i= j-1; i > 0; i-- ) |
---|
653 | *a[i] += (*a[i+1]*x); |
---|
654 | for (i= 0; i < j; i++ ) |
---|
655 | *a[i] = *a[i+1]; |
---|
656 | } |
---|
657 | else |
---|
658 | { |
---|
659 | gmp_complex y(o/x); |
---|
660 | for (i= 1; i < j; i++) |
---|
661 | *a[i] += (*a[i-1]*y); |
---|
662 | } |
---|
663 | } |
---|
664 | |
---|
665 | void rootContainer::divquad(gmp_complex **a, gmp_complex x, int j) |
---|
666 | { |
---|
667 | int i; |
---|
668 | gmp_float o(1.0),p(x.real()+x.real()), |
---|
669 | q((x.real()*x.real())+(x.imag()*x.imag())); |
---|
670 | |
---|
671 | if (abs(x)<o) |
---|
672 | { |
---|
673 | *a[j-1] += (*a[j]*p); |
---|
674 | for (i= j-2; i > 1; i-- ) |
---|
675 | *a[i] += ((*a[i+1]*p)-(*a[i+2]*q)); |
---|
676 | for (i= 0; i < j-1; i++ ) |
---|
677 | *a[i] = *a[i+2]; |
---|
678 | } |
---|
679 | else |
---|
680 | { |
---|
681 | p = p/q; |
---|
682 | q = o/q; |
---|
683 | *a[1] += (*a[0]*p); |
---|
684 | for (i= 2; i < j-1; i++) |
---|
685 | *a[i] += ((*a[i-1]*p)-(*a[i-2]*q)); |
---|
686 | } |
---|
687 | } |
---|
688 | |
---|
689 | void rootContainer::solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j) |
---|
690 | { |
---|
691 | gmp_float zero(0.0); |
---|
692 | |
---|
693 | if ((j>k) |
---|
694 | &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero()))) |
---|
695 | { |
---|
696 | gmp_complex sq(zero); |
---|
697 | gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]); |
---|
698 | gmp_complex disk((h1 * h1) - h2); |
---|
699 | if (disk.imag()==zero) |
---|
700 | { |
---|
701 | if (disk.real()<zero) |
---|
702 | { |
---|
703 | sq.real(zero); |
---|
704 | sq.imag(sqrt(-disk.real())); |
---|
705 | } |
---|
706 | else |
---|
707 | sq = (gmp_complex)sqrt(disk.real()); |
---|
708 | } |
---|
709 | else |
---|
710 | sq = sqrt(disk); |
---|
711 | *r[k+1] = sq - h1; |
---|
712 | sq += h1; |
---|
713 | *r[k] = (gmp_complex)0.0-sq; |
---|
714 | if(sq.imag()==zero) |
---|
715 | { |
---|
716 | k = j; |
---|
717 | j++; |
---|
718 | } |
---|
719 | else |
---|
720 | { |
---|
721 | j = k; |
---|
722 | k--; |
---|
723 | } |
---|
724 | } |
---|
725 | else |
---|
726 | { |
---|
727 | if (((*a[1]).real()==zero) && ((*a[1]).imag()==zero)) |
---|
728 | { |
---|
729 | WerrorS("precision lost, try again with higher precision"); |
---|
730 | } |
---|
731 | else |
---|
732 | { |
---|
733 | *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]); |
---|
734 | if(r[k]->imag()==zero) |
---|
735 | j++; |
---|
736 | else |
---|
737 | k--; |
---|
738 | } |
---|
739 | } |
---|
740 | } |
---|
741 | |
---|
742 | void rootContainer::sortroots(gmp_complex **ro, int r, int c, bool isf) |
---|
743 | { |
---|
744 | int j; |
---|
745 | |
---|
746 | for (j=0; j<r; j++) // the real roots |
---|
747 | sortre(ro, j, r, 1); |
---|
748 | if (c>=tdg) return; |
---|
749 | if (isf) |
---|
750 | { |
---|
751 | for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly |
---|
752 | sortre(ro, j, tdg-1, 2); |
---|
753 | } |
---|
754 | else |
---|
755 | { |
---|
756 | for (j=c; j+1<tdg; j++) // the complex roots for a general poly |
---|
757 | sortre(ro, j, tdg-1, 1); |
---|
758 | } |
---|
759 | } |
---|
760 | |
---|
761 | void rootContainer::sortre(gmp_complex **r, int l, int u, int inc) |
---|
762 | { |
---|
763 | int pos,i; |
---|
764 | gmp_complex *x,*y; |
---|
765 | |
---|
766 | pos = l; |
---|
767 | x = r[pos]; |
---|
768 | for (i=l+inc; i<=u; i+=inc) |
---|
769 | { |
---|
770 | if (r[i]->real()<x->real()) |
---|
771 | { |
---|
772 | pos = i; |
---|
773 | x = r[pos]; |
---|
774 | } |
---|
775 | } |
---|
776 | if (pos>l) |
---|
777 | { |
---|
778 | if (inc==1) |
---|
779 | { |
---|
780 | for (i=pos; i>l; i--) |
---|
781 | r[i] = r[i-1]; |
---|
782 | r[l] = x; |
---|
783 | } |
---|
784 | else |
---|
785 | { |
---|
786 | y = r[pos+1]; |
---|
787 | for (i=pos+1; i+1>l; i--) |
---|
788 | r[i] = r[i-2]; |
---|
789 | if (x->imag()>y->imag()) |
---|
790 | { |
---|
791 | r[l] = x; |
---|
792 | r[l+1] = y; |
---|
793 | } |
---|
794 | else |
---|
795 | { |
---|
796 | r[l] = y; |
---|
797 | r[l+1] = x; |
---|
798 | } |
---|
799 | } |
---|
800 | } |
---|
801 | else if ((inc==2)&&(x->imag()<r[l+1]->imag())) |
---|
802 | { |
---|
803 | r[l] = r[l+1]; |
---|
804 | r[l+1] = x; |
---|
805 | } |
---|
806 | } |
---|
807 | |
---|
808 | void rootContainer::computefx(gmp_complex **a, gmp_complex x, int m, |
---|
809 | gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, |
---|
810 | gmp_float &ex, gmp_float &ef) |
---|
811 | { |
---|
812 | int k; |
---|
813 | |
---|
814 | f0= *a[m]; |
---|
815 | ef= abs(f0); |
---|
816 | f1= gmp_complex(0.0); |
---|
817 | f2= f1; |
---|
818 | ex= abs(x); |
---|
819 | |
---|
820 | for ( k= m-1; k >= 0; k-- ) |
---|
821 | { |
---|
822 | f2 = ( x * f2 ) + f1; |
---|
823 | f1 = ( x * f1 ) + f0; |
---|
824 | f0 = ( x * f0 ) + *a[k]; |
---|
825 | ef = abs( f0 ) + ( ex * ef ); |
---|
826 | } |
---|
827 | } |
---|
828 | |
---|
829 | void rootContainer::computegx(gmp_complex **a, gmp_complex x, int m, |
---|
830 | gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, |
---|
831 | gmp_float &ex, gmp_float &ef) |
---|
832 | { |
---|
833 | int k; |
---|
834 | |
---|
835 | f0= *a[0]; |
---|
836 | ef= abs(f0); |
---|
837 | f1= gmp_complex(0.0); |
---|
838 | f2= f1; |
---|
839 | ex= abs(x); |
---|
840 | |
---|
841 | for ( k= 1; k <= m; k++ ) |
---|
842 | { |
---|
843 | f2 = ( x * f2 ) + f1; |
---|
844 | f1 = ( x * f1 ) + f0; |
---|
845 | f0 = ( x * f0 ) + *a[k]; |
---|
846 | ef = abs( f0 ) + ( ex * ef ); |
---|
847 | } |
---|
848 | } |
---|
849 | |
---|
850 | //----------------------------------------------------------------------------- |
---|
851 | //-------------- rootArranger ------------------------------------------------- |
---|
852 | //----------------------------------------------------------------------------- |
---|
853 | |
---|
854 | //-> rootArranger::rootArranger(...) |
---|
855 | rootArranger::rootArranger( rootContainer ** _roots, |
---|
856 | rootContainer ** _mu, |
---|
857 | const int _howclean ) |
---|
858 | : roots(_roots), mu(_mu), howclean(_howclean) |
---|
859 | { |
---|
860 | found_roots=false; |
---|
861 | } |
---|
862 | //<- |
---|
863 | |
---|
864 | //-> void rootArranger::solve_all() |
---|
865 | void rootArranger::solve_all() |
---|
866 | { |
---|
867 | int i; |
---|
868 | found_roots= true; |
---|
869 | |
---|
870 | // find roots of polys given by coeffs in roots |
---|
871 | rc= roots[0]->getAnzElems(); |
---|
872 | for ( i= 0; i < rc; i++ ) |
---|
873 | if ( !roots[i]->solver( howclean ) ) |
---|
874 | { |
---|
875 | found_roots= false; |
---|
876 | return; |
---|
877 | } |
---|
878 | // find roots of polys given by coeffs in mu |
---|
879 | mc= mu[0]->getAnzElems(); |
---|
880 | for ( i= 0; i < mc; i++ ) |
---|
881 | if ( ! mu[i]->solver( howclean ) ) |
---|
882 | { |
---|
883 | found_roots= false; |
---|
884 | return; |
---|
885 | } |
---|
886 | } |
---|
887 | //<- |
---|
888 | |
---|
889 | //-> void rootArranger::arrange() |
---|
890 | void rootArranger::arrange() |
---|
891 | { |
---|
892 | gmp_complex tmp,zwerg; |
---|
893 | int anzm= mu[0]->getAnzElems(); |
---|
894 | int anzr= roots[0]->getAnzRoots(); |
---|
895 | int xkoord, r, rtest, xk, mtest; |
---|
896 | bool found; |
---|
897 | //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5)); |
---|
898 | gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3))); |
---|
899 | |
---|
900 | for ( xkoord= 0; xkoord < anzm; xkoord++ ) { // für x1,x2, x1,x2,x3, x1,x2,...,xn |
---|
901 | for ( r= 0; r < anzr; r++ ) { // für jede Nullstelle |
---|
902 | // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] + |
---|
903 | // ... + (xkoord-koordinate) * evp[xkoord] |
---|
904 | tmp= gmp_complex(); |
---|
905 | for ( xk =0; xk <= xkoord; xk++ ) |
---|
906 | { |
---|
907 | tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1 |
---|
908 | } |
---|
909 | found= false; |
---|
910 | for ( rtest= r; rtest < anzr; rtest++ ) { // für jede Nullstelle |
---|
911 | zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2 |
---|
912 | for ( mtest= 0; mtest < anzr; mtest++ ) |
---|
913 | { |
---|
914 | // if ( tmp == (*mu[xkoord])[mtest] ) |
---|
915 | // { |
---|
916 | if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) && |
---|
917 | (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) && |
---|
918 | ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) && |
---|
919 | (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) ) |
---|
920 | { |
---|
921 | roots[xk]->swapRoots( r, rtest ); |
---|
922 | found= true; |
---|
923 | break; |
---|
924 | } |
---|
925 | } |
---|
926 | } // rtest |
---|
927 | if ( !found ) |
---|
928 | { |
---|
929 | Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r); |
---|
930 | #ifdef mprDEBUG_PROT |
---|
931 | WarnS("One of these ..."); |
---|
932 | for ( rtest= r; rtest < anzr; rtest++ ) |
---|
933 | { |
---|
934 | tmp= gmp_complex(); |
---|
935 | for ( xk =0; xk <= xkoord; xk++ ) |
---|
936 | { |
---|
937 | tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); |
---|
938 | } |
---|
939 | tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2 |
---|
940 | Warn(" %s",complexToStr(tmp,gmp_output_digits+1),rtest); |
---|
941 | } |
---|
942 | WarnS(" ... must match to one of these:"); |
---|
943 | for ( mtest= 0; mtest < anzr; mtest++ ) |
---|
944 | { |
---|
945 | Warn(" %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1)); |
---|
946 | } |
---|
947 | #endif |
---|
948 | } |
---|
949 | } // r |
---|
950 | } // xkoord |
---|
951 | } |
---|
952 | //<- |
---|
953 | |
---|
954 | //----------------------------------------------------------------------------- |
---|
955 | //-------------- simplex ----- ------------------------------------------------ |
---|
956 | //----------------------------------------------------------------------------- |
---|
957 | |
---|
958 | // #ifdef mprDEBUG_PROT |
---|
959 | // #define error(a) a |
---|
960 | // #else |
---|
961 | // #define error(a) |
---|
962 | // #endif |
---|
963 | |
---|
964 | #define error(a) a |
---|
965 | |
---|
966 | #define MAXPOINTS 1000 |
---|
967 | |
---|
968 | //-> simplex::* |
---|
969 | // |
---|
970 | simplex::simplex( int rows, int cols ) |
---|
971 | : LiPM_cols(cols), LiPM_rows(rows) |
---|
972 | { |
---|
973 | int i; |
---|
974 | |
---|
975 | LiPM_rows=LiPM_rows+3; |
---|
976 | LiPM_cols=LiPM_cols+2; |
---|
977 | |
---|
978 | LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) ); // LP matrix |
---|
979 | for( i= 0; i < LiPM_rows; i++ ) |
---|
980 | { |
---|
981 | // Mem must be allocated aligned, also for type double! |
---|
982 | LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) ); |
---|
983 | } |
---|
984 | |
---|
985 | iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) ); |
---|
986 | izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) ); |
---|
987 | |
---|
988 | m=n=m1=m2=m3=icase=0; |
---|
989 | |
---|
990 | #ifdef mprDEBUG_ALL |
---|
991 | Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols); |
---|
992 | #endif |
---|
993 | } |
---|
994 | |
---|
995 | simplex::~simplex() |
---|
996 | { |
---|
997 | // clean up |
---|
998 | int i; |
---|
999 | for( i= 0; i < LiPM_rows; i++ ) |
---|
1000 | { |
---|
1001 | omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) ); |
---|
1002 | } |
---|
1003 | omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) ); |
---|
1004 | |
---|
1005 | omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) ); |
---|
1006 | omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) ); |
---|
1007 | } |
---|
1008 | |
---|
1009 | BOOLEAN simplex::mapFromMatrix( matrix m ) |
---|
1010 | { |
---|
1011 | int i,j; |
---|
1012 | // if ( MATROWS( m ) > LiPM_rows || MATCOLS( m ) > LiPM_cols ) { |
---|
1013 | // WarnS(""); |
---|
1014 | // return FALSE; |
---|
1015 | // } |
---|
1016 | |
---|
1017 | number coef; |
---|
1018 | for ( i= 1; i <= MATROWS( m ); i++ ) |
---|
1019 | { |
---|
1020 | for ( j= 1; j <= MATCOLS( m ); j++ ) |
---|
1021 | { |
---|
1022 | if ( MATELEM(m,i,j) != NULL ) |
---|
1023 | { |
---|
1024 | coef= pGetCoeff( MATELEM(m,i,j) ); |
---|
1025 | if ( coef != NULL && !nIsZero(coef) ) |
---|
1026 | LiPM[i][j]= (double)(*(gmp_float*)coef); |
---|
1027 | //#ifdef mpr_DEBUG_PROT |
---|
1028 | //Print("%f ",LiPM[i][j]); |
---|
1029 | //#endif |
---|
1030 | } |
---|
1031 | } |
---|
1032 | // PrintLn(); |
---|
1033 | } |
---|
1034 | |
---|
1035 | return TRUE; |
---|
1036 | } |
---|
1037 | |
---|
1038 | matrix simplex::mapToMatrix( matrix m ) |
---|
1039 | { |
---|
1040 | int i,j; |
---|
1041 | // if ( MATROWS( m ) < LiPM_rows-3 || MATCOLS( m ) < LiPM_cols-2 ) { |
---|
1042 | // WarnS(""); |
---|
1043 | // return NULL; |
---|
1044 | // } |
---|
1045 | |
---|
1046 | //Print(" %d x %d\n",MATROWS( m ),MATCOLS( m )); |
---|
1047 | |
---|
1048 | number coef; |
---|
1049 | gmp_float * bla; |
---|
1050 | for ( i= 1; i <= MATROWS( m ); i++ ) |
---|
1051 | { |
---|
1052 | for ( j= 1; j <= MATCOLS( m ); j++ ) |
---|
1053 | { |
---|
1054 | pDelete( &(MATELEM(m,i,j)) ); |
---|
1055 | MATELEM(m,i,j)= NULL; |
---|
1056 | //Print(" %3.0f ",LiPM[i][j]); |
---|
1057 | if ( LiPM[i][j] != 0.0 ) |
---|
1058 | { |
---|
1059 | bla= new gmp_float(LiPM[i][j]); |
---|
1060 | coef= (number)bla; |
---|
1061 | MATELEM(m,i,j)= pOne(); |
---|
1062 | pSetCoeff( MATELEM(m,i,j), coef ); |
---|
1063 | } |
---|
1064 | } |
---|
1065 | //PrintLn(); |
---|
1066 | } |
---|
1067 | |
---|
1068 | return m; |
---|
1069 | } |
---|
1070 | |
---|
1071 | intvec * simplex::posvToIV() |
---|
1072 | { |
---|
1073 | int i; |
---|
1074 | intvec * iv = new intvec( m ); |
---|
1075 | for ( i= 1; i <= m; i++ ) |
---|
1076 | { |
---|
1077 | IMATELEM(*iv,i,1)= iposv[i]; |
---|
1078 | } |
---|
1079 | return iv; |
---|
1080 | } |
---|
1081 | |
---|
1082 | intvec * simplex::zrovToIV() |
---|
1083 | { |
---|
1084 | int i; |
---|
1085 | intvec * iv = new intvec( n ); |
---|
1086 | for ( i= 1; i <= n; i++ ) |
---|
1087 | { |
---|
1088 | IMATELEM(*iv,i,1)= izrov[i]; |
---|
1089 | } |
---|
1090 | return iv; |
---|
1091 | } |
---|
1092 | |
---|
1093 | void simplex::compute() |
---|
1094 | { |
---|
1095 | int i,ip,ir,is,k,kh,kp,m12,nl1,nl2; |
---|
1096 | int *l1,*l2,*l3; |
---|
1097 | mprfloat q1, bmax; |
---|
1098 | |
---|
1099 | if ( m != (m1+m2+m3) ) |
---|
1100 | { |
---|
1101 | // error: bad input |
---|
1102 | error(WarnS("simplex::compute: Bad input constraint counts!");) |
---|
1103 | icase=-2; |
---|
1104 | return; |
---|
1105 | } |
---|
1106 | |
---|
1107 | l1= (int *) omAlloc0( (n+1) * sizeof(int) ); |
---|
1108 | l2= (int *) omAlloc0( (m+1) * sizeof(int) ); |
---|
1109 | l3= (int *) omAlloc0( (m+1) * sizeof(int) ); |
---|
1110 | |
---|
1111 | nl1= n; |
---|
1112 | for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k; |
---|
1113 | nl2=m; |
---|
1114 | for ( i=1; i<=m; i++ ) |
---|
1115 | { |
---|
1116 | if ( LiPM[i+1][1] < 0.0 ) |
---|
1117 | { |
---|
1118 | // error: bad input |
---|
1119 | error(WarnS("simplex::compute: Bad input tableau!");) |
---|
1120 | error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);) |
---|
1121 | icase=-2; |
---|
1122 | // free mem l1,l2,l3; |
---|
1123 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
1124 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
1125 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
1126 | return; |
---|
1127 | } |
---|
1128 | l2[i]= i; |
---|
1129 | iposv[i]= n+i; |
---|
1130 | } |
---|
1131 | for ( i=1; i<=m2; i++) l3[i]= 1; |
---|
1132 | ir= 0; |
---|
1133 | if (m2+m3) |
---|
1134 | { |
---|
1135 | ir=1; |
---|
1136 | for ( k=1; k <= (n+1); k++ ) |
---|
1137 | { |
---|
1138 | q1=0.0; |
---|
1139 | for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k]; |
---|
1140 | LiPM[m+2][k]= -q1; |
---|
1141 | } |
---|
1142 | |
---|
1143 | do |
---|
1144 | { |
---|
1145 | simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax); |
---|
1146 | if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS ) |
---|
1147 | { |
---|
1148 | icase= -1; // no solution found |
---|
1149 | // free mem l1,l2,l3; |
---|
1150 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
1151 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
1152 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
1153 | return; |
---|
1154 | } |
---|
1155 | else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS ) |
---|
1156 | { |
---|
1157 | m12= m1+m2+1; |
---|
1158 | if ( m12 <= m ) |
---|
1159 | { |
---|
1160 | for ( ip= m12; ip <= m; ip++ ) |
---|
1161 | { |
---|
1162 | if ( iposv[ip] == (ip+n) ) |
---|
1163 | { |
---|
1164 | simp1(LiPM,ip,l1,nl1,1,&kp,&bmax); |
---|
1165 | if ( fabs(bmax) >= SIMPLEX_EPS) |
---|
1166 | goto one; |
---|
1167 | } |
---|
1168 | } |
---|
1169 | } |
---|
1170 | ir= 0; |
---|
1171 | --m12; |
---|
1172 | if ( m1+1 <= m12 ) |
---|
1173 | for ( i=m1+1; i <= m12; i++ ) |
---|
1174 | if ( l3[i-m1] == 1 ) |
---|
1175 | for ( k=1; k <= n+1; k++ ) |
---|
1176 | LiPM[i+1][k] = -(LiPM[i+1][k]); |
---|
1177 | break; |
---|
1178 | } |
---|
1179 | //#if DEBUG |
---|
1180 | //print_bmat( a, m+2, n+3); |
---|
1181 | //#endif |
---|
1182 | simp2(LiPM,n,l2,nl2,&ip,kp,&q1); |
---|
1183 | if ( ip == 0 ) |
---|
1184 | { |
---|
1185 | icase = -1; // no solution found |
---|
1186 | // free mem l1,l2,l3; |
---|
1187 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
1188 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
1189 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
1190 | return; |
---|
1191 | } |
---|
1192 | one: simp3(LiPM,m+1,n,ip,kp); |
---|
1193 | // #if DEBUG |
---|
1194 | // print_bmat(a,m+2,n+3); |
---|
1195 | // #endif |
---|
1196 | if ( iposv[ip] >= (n+m1+m2+1)) |
---|
1197 | { |
---|
1198 | for ( k= 1; k <= nl1; k++ ) |
---|
1199 | if ( l1[k] == kp ) break; |
---|
1200 | --nl1; |
---|
1201 | for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1]; |
---|
1202 | ++(LiPM[m+2][kp+1]); |
---|
1203 | for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]); |
---|
1204 | } |
---|
1205 | else |
---|
1206 | { |
---|
1207 | if ( iposv[ip] >= (n+m1+1) ) |
---|
1208 | { |
---|
1209 | kh= iposv[ip]-m1-n; |
---|
1210 | if ( l3[kh] ) |
---|
1211 | { |
---|
1212 | l3[kh]= 0; |
---|
1213 | ++(LiPM[m+2][kp+1]); |
---|
1214 | for ( i=1; i<= m+2; i++ ) |
---|
1215 | LiPM[i][kp+1] = -(LiPM[i][kp+1]); |
---|
1216 | } |
---|
1217 | } |
---|
1218 | } |
---|
1219 | is= izrov[kp]; |
---|
1220 | izrov[kp]= iposv[ip]; |
---|
1221 | iposv[ip]= is; |
---|
1222 | } while (ir); |
---|
1223 | } |
---|
1224 | /* end of phase 1, have feasible sol, now optimize it */ |
---|
1225 | loop |
---|
1226 | { |
---|
1227 | // #if DEBUG |
---|
1228 | // print_bmat( a, m+1, n+5); |
---|
1229 | // #endif |
---|
1230 | simp1(LiPM,0,l1,nl1,0,&kp,&bmax); |
---|
1231 | if (bmax <= /*SIMPLEX_EPS*/0.0) |
---|
1232 | { |
---|
1233 | icase=0; // finite solution found |
---|
1234 | // free mem l1,l2,l3 |
---|
1235 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
1236 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
1237 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
1238 | return; |
---|
1239 | } |
---|
1240 | simp2(LiPM,n,l2,nl2,&ip,kp,&q1); |
---|
1241 | if (ip == 0) |
---|
1242 | { |
---|
1243 | //printf("Unbounded:"); |
---|
1244 | // #if DEBUG |
---|
1245 | // print_bmat( a, m+1, n+1); |
---|
1246 | // #endif |
---|
1247 | icase=1; /* unbounded */ |
---|
1248 | // free mem |
---|
1249 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
1250 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
1251 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
1252 | return; |
---|
1253 | } |
---|
1254 | simp3(LiPM,m,n,ip,kp); |
---|
1255 | is= izrov[kp]; |
---|
1256 | izrov[kp]= iposv[ip]; |
---|
1257 | iposv[ip]= is; |
---|
1258 | }/*for ;;*/ |
---|
1259 | } |
---|
1260 | |
---|
1261 | void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax ) |
---|
1262 | { |
---|
1263 | int k; |
---|
1264 | mprfloat test; |
---|
1265 | |
---|
1266 | if( nll <= 0) |
---|
1267 | { /* init'tion: fixed */ |
---|
1268 | *bmax = 0.0; |
---|
1269 | return; |
---|
1270 | } |
---|
1271 | *kp=ll[1]; |
---|
1272 | *bmax=a[mm+1][*kp+1]; |
---|
1273 | for (k=2;k<=nll;k++) |
---|
1274 | { |
---|
1275 | if (iabf == 0) |
---|
1276 | { |
---|
1277 | test=a[mm+1][ll[k]+1]-(*bmax); |
---|
1278 | if (test > 0.0) |
---|
1279 | { |
---|
1280 | *bmax=a[mm+1][ll[k]+1]; |
---|
1281 | *kp=ll[k]; |
---|
1282 | } |
---|
1283 | } |
---|
1284 | else |
---|
1285 | { /* abs values: have fixed it */ |
---|
1286 | test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax); |
---|
1287 | if (test > 0.0) |
---|
1288 | { |
---|
1289 | *bmax=a[mm+1][ll[k]+1]; |
---|
1290 | *kp=ll[k]; |
---|
1291 | } |
---|
1292 | } |
---|
1293 | } |
---|
1294 | } |
---|
1295 | |
---|
1296 | void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 ) |
---|
1297 | { |
---|
1298 | int k,ii,i; |
---|
1299 | mprfloat qp,q0,q; |
---|
1300 | |
---|
1301 | *ip= 0; |
---|
1302 | for ( i=1; i <= nl2; i++ ) |
---|
1303 | { |
---|
1304 | if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS ) |
---|
1305 | { |
---|
1306 | *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1]; |
---|
1307 | *ip= l2[i]; |
---|
1308 | for ( i= i+1; i <= nl2; i++ ) |
---|
1309 | { |
---|
1310 | ii= l2[i]; |
---|
1311 | if (a[ii+1][kp+1] < -SIMPLEX_EPS) |
---|
1312 | { |
---|
1313 | q= -a[ii+1][1] / a[ii+1][kp+1]; |
---|
1314 | if (q - *q1 < -SIMPLEX_EPS) |
---|
1315 | { |
---|
1316 | *ip=ii; |
---|
1317 | *q1=q; |
---|
1318 | } |
---|
1319 | else if (q - *q1 < SIMPLEX_EPS) |
---|
1320 | { |
---|
1321 | for ( k=1; k<= n; k++ ) |
---|
1322 | { |
---|
1323 | qp= -a[*ip+1][k+1]/a[*ip+1][kp+1]; |
---|
1324 | q0= -a[ii+1][k+1]/a[ii+1][kp+1]; |
---|
1325 | if ( q0 != qp ) break; |
---|
1326 | } |
---|
1327 | if ( q0 < qp ) *ip= ii; |
---|
1328 | } |
---|
1329 | } |
---|
1330 | } |
---|
1331 | } |
---|
1332 | } |
---|
1333 | } |
---|
1334 | |
---|
1335 | void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp ) |
---|
1336 | { |
---|
1337 | int kk,ii; |
---|
1338 | mprfloat piv; |
---|
1339 | |
---|
1340 | piv= 1.0 / a[ip+1][kp+1]; |
---|
1341 | for ( ii=1; ii <= i1+1; ii++ ) |
---|
1342 | { |
---|
1343 | if ( ii -1 != ip ) |
---|
1344 | { |
---|
1345 | a[ii][kp+1] *= piv; |
---|
1346 | for ( kk=1; kk <= k1+1; kk++ ) |
---|
1347 | if ( kk-1 != kp ) |
---|
1348 | a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1]; |
---|
1349 | } |
---|
1350 | } |
---|
1351 | for ( kk=1; kk<= k1+1; kk++ ) |
---|
1352 | if ( kk-1 != kp ) a[ip+1][kk] *= -piv; |
---|
1353 | a[ip+1][kp+1]= piv; |
---|
1354 | } |
---|
1355 | //<- |
---|
1356 | |
---|
1357 | //----------------------------------------------------------------------------- |
---|
1358 | |
---|
1359 | //#endif // HAVE_MPR |
---|
1360 | |
---|
1361 | // local Variables: *** |
---|
1362 | // folded-file: t *** |
---|
1363 | // compile-command-1: "make installg" *** |
---|
1364 | // compile-command-2: "make install" *** |
---|
1365 | // End: *** |
---|